Search code examples
rlme4anovaxtable

How to convert Afex or car ANOVA models to lmer? Observed variables


In the afex package we can find this example of ANOVA analysis:

data(obk.long, package = "afex")
# estimate mixed ANOVA on the full design:
# can be written in any of these ways: 

aov_car(value ~ treatment * gender + Error(id/(phase*hour)), data = obk.long,
        observed = "gender")

aov_4(value ~ treatment * gender + (phase*hour|id), data = obk.long,
        observed = "gender")

aov_ez("id", "value", obk.long, between = c("treatment", "gender"), 
       within = c("phase", "hour"), observed = "gender")

My question is, How can I write the same model in lme4? In particular, I don't know how to include the "observed" term?

If I just write

lmer(value ~ treatment * gender + (phase*hour|id), data = obk.long,
     observed = "gender")

I get an error telling that observed is not a valid option.

Furthermore, if I just remove the observed option lmer produces the error:

Error: number of observations (=240) <= number of random effects (=240) for term (phase * hour | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable.

Where in the lmer syntax do I specify the "between" or "within" variable?. As far as I know you just write the dependent variable on the left side and all other variables on the right side, and the error term as (1|id).

The package "car" uses the idata for the intra-subject variable.


Solution

  • I might not know enough about classical ANOVA theory to answer this question completely, but I'll take a crack. First, a couple of points:

    • the observed argument appears only to be relevant for the computation of effect size.

    observed: ‘character’ vector indicating which of the variables are observed (i.e, measured) as compared to experimentally manipulated. The default effect size reported (generalized eta-squared) requires correct specification of the obsered [sic] (in contrast to manipulated) variables.

    ... so I think you'd be safe leaving it out.

    • if you want to override the error you can use
      control=lmerControl(check.nobs.vs.nRE="ignore")
    

    ... but this probably isn't the right way forward.

    I think but am not sure that this is the right way:

    m1 <- lmer(value ~ treatment * gender + (1|id/phase:hour), data = obk.long,
        control=lmerControl(check.nobs.vs.nRE="ignore",
                                check.nobs.vs.nlev="ignore"),
        contrasts=list(treatment=contr.sum,gender=contr.sum))
    

    This specifies that the interaction of phase and hour varies within id. The residual variance and (phase by hour within id) variance are confounded (which is why we need the overriding lmerControl() specification), so don't trust those particular variance estimates. However, the main effects of treatment and gender should be handled just the same. If you load lmerTest instead of lmer and run summary(m1) or anova(m1) it gives you the same degrees of freedom (10) for the fixed (gender and treatment) effects that are computed by afex.

    lme gives comparable answers, but needs to have the phase-by-hour interaction constructed beforehand:

    library(nlme)
    obk.long$ph <- with(obk.long,interaction(phase,hour))
    m2 <- lme(value ~ treatment * gender,
                 random=~1|id/ph, data = obk.long,
        contrasts=list(treatment=contr.sum,gender=contr.sum))
    anova(m2,type="marginal")
    

    I don't know how to reconstruct afex's tests of the random effects.