CODEWITHSUNDEEP
rlme4anovaxtable
skan
skan

Reputation: 7720

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.

Upvotes: 6

Views: 952

Answers (2)

Henrik
Henrik

Reputation: 14450

As Ben Bolker correctly says, simply leave observed out.

Furthermore, I would not recommend to do what you want to do. Using a mixed model for a data set without replications within each cell of the design per participant is somewhat questionable as it is not really clear how to specify the random effects structure. Importantly, the Barr et al. maxim of "keep it maximal" does not work here as you realized. The problem is that the model is overparametrized (hence the error from lmer).

I recommend using the ANOVA. More discussion on exactly this question can be found on a crossvalidated thread where Ben and me discussed this more thoroughly.

Upvotes: 4

Ben Bolker
Ben Bolker

Reputation: 226162

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.

Upvotes: 4

Related Questions