AnovaMixedModels

The implementation of ANOVA for mixed-effects models is primarily based on MixedModels. The syntax is similar to anova for GLM.

using AnovaMixedModels

Linear mixed-effects model

We get a dataset from R directly by RCall.

R"""data("anxiety", package = "datarium")"""
anxiety = stack(rcopy(R"anxiety"), [:t1, :t2, :t3], [:id, :group], variable_name = :time, value_name = :score)
anxiety = combine(anxiety, Not(:time), :time => ByRow(x->parse(Int, replace(String(x), "t"=>""))) => :time)

We can fit a linear mixed-effects model first. lme is an alias for fit(LinearMixedModel, formula, data, args...).

lmm1 = lme(@formula(score ~ group * time + (1|id)), anxiety)
anova(lmm1)

Analysis of Variance

Type 1 test / F test

score ~ 1 + group + time + group & time + (1 | id)

Table:
───────────────────────────────────────────────
              DOF  Res.DOF    F value  Pr(>|F|)
───────────────────────────────────────────────
(Intercept)     1       87  5019.3886    <1e-77
group           2       42     4.4554    0.0176
time            1       87   604.7798    <1e-40
group & time    2       87   159.3211    <1e-29
───────────────────────────────────────────────

Alternatively, we can use anova_lme. Like anova_lm, this function will fit and store a model; in this case, a LinearMixedModel fitted by Restricted maximum likelihood.

aov = anova_lme(@formula(score ~ group * time + (1|id)), anxiety, type = 3)

Analysis of Variance

Type 3 test / F test

score ~ 1 + group + time + group & time + (1 | id)

Table:
───────────────────────────────────────────────
              DOF  Res.DOF    F value  Pr(>|F|)
───────────────────────────────────────────────
(Intercept)     1       87  1756.6506    <1e-58
group           2       42     3.1340    0.0539
time            1       87    23.2498    <1e-05
group & time    2       87   161.1736    <1e-29
───────────────────────────────────────────────

aov.anovamodel.model.optsum.REML

true

For likeihood-ratio test, all submodels are fitted. The model should be fitted by maximum likelihood estimation.

anova(LRT, lmm1)

Analysis of Variance

Type 1 test / Likelihood-ratio test

Model 1: score ~ 0 + (1 | id)
Model 2: score ~ 1 + (1 | id)
Model 3: score ~ 1 + group + (1 | id)
Model 4: score ~ 1 + group + time + (1 | id)
Model 5: score ~ 1 + group + time + group & time + (1 | id)

Table:
────────────────────────────────────────────────────
   DOF  ΔDOF  Res.DOF  Deviance        χ²  Pr(>|χ²|)
────────────────────────────────────────────────────
1    2            133    701.74               
2    3     1      132    495.55  206.1822     <1e-46
3    5     2      130    487.08    8.4747     0.0144
4    6     1      129    404.81   82.2717     <1e-18
5    8     2      127    265.43  139.3790     <1e-30
────────────────────────────────────────────────────

When comparing multiple mixed models, likelihood-ratio test is used by default. It's also identical to StatsModels.lrtest and MixedModels.likelihoodratiotest.

lmms = nestedmodels(lmm1)
anova(lmms) # == anova(LRT, lmm1)

Analysis of Variance

Type 1 test / Likelihood-ratio test

Model 1: score ~ 0 + (1 | id)
Model 2: score ~ 1 + (1 | id)
Model 3: score ~ 1 + group + (1 | id)
Model 4: score ~ 1 + group + time + (1 | id)
Model 5: score ~ 1 + group + time + group & time + (1 | id)

Table:
────────────────────────────────────────────────────
   DOF  ΔDOF  Res.DOF  Deviance        χ²  Pr(>|χ²|)
────────────────────────────────────────────────────
1    2            133    701.74               
2    3     1      132    495.55  206.1822     <1e-46
3    5     2      130    487.08    8.4747     0.0144
4    6     1      129    404.81   82.2717     <1e-18
5    8     2      127    265.43  139.3790     <1e-30
────────────────────────────────────────────────────

MixedModels.likelihoodratiotest(lmms.model[2:end]...)

Comparing between LinearModel and LinearMixedModel is also available.

lm1 = lm(@formula(score ~ group * time), anxiety)
lmm2 = lme(@formula(score ~ group * time + (group|id)), anxiety)
anova(lm1, lmm1, lmm2)

Analysis of Variance

Type 1 test / Likelihood-ratio test

Model 1: score ~ 1 + group + time + group & time
Model 2: score ~ 1 + group + time + group & time + (1 | id)
Model 3: score ~ 1 + group + time + group & time + (1 + group | id)

Table:
─────────────────────────────────────────────────────────────────────
   DOF  ΔDOF  Res.DOF     AIC     BIC  -2 logLik        χ²  Pr(>|χ²|)
─────────────────────────────────────────────────────────────────────
1    7            129  508.72  529.06     494.72               
2    8     1      127  281.43  304.67     265.43  229.2935     <1e-51
3   13     5      122  290.79  328.56     264.79    0.6349     0.9864
─────────────────────────────────────────────────────────────────────

Generalized linear mixed-effects model

The following is an example of generalized mixed model. glme is an alias for fit(GeneralizedLinearMixedModel, formula, data, args...).

R"""data("toenail", package = "HSAUR2")"""
toenail = rcopy(R"toenail")

glmm1 = glme(@formula(outcome ~ visit + treatment + (1|patientID)), toenail, Binomial(), LogitLink(), nAGQ = 20, wts = ones(Float64, size(toenail, 1)));
glmm2 = glme(@formula(outcome ~ visit * treatment + (1|patientID)), toenail, Binomial(), LogitLink(), nAGQ = 20, wts = ones(Float64, size(toenail, 1)));
anova(glmm1, glmm2)

Analysis of Variance

Type 1 test / Likelihood-ratio test

Model 1: outcome ~ 1 + visit + treatment + (1 | patientID)
Model 2: outcome ~ 1 + visit + treatment + visit & treatment + (1 | patientID)

Table:
──────────────────────────────────────────────────
   DOF  ΔDOF  Res.DOF  Deviance      χ²  Pr(>|χ²|)
──────────────────────────────────────────────────
1    4           1904   1246.12             
2    5     1     1903   1242.39  3.7270     0.0535
──────────────────────────────────────────────────