R: Mixed ANOVA works with the aov{stats} function and fails with the anova_test{rstatix} function. Why?

Question

I can implement a rather complicated mixed ANOVA in R with the aov{stats} or the lmer{lme4} functions, but not with the anova_test{rstatix} nor the Anova{car} functions. My question follows: How can I make the anova_test{rstatix} or Anova{car} functions work to access Greenhouse-Geisser corrections?

I have a factorial, incomplete blocks experiment with repeated measures. I replicated factorial nitrogen (N), phosphorus (P) and potassium (K) addition four times, blocking control, +NP, +NK and +PK treatments versus +N, +P, +K and +NPK treatments within each replicate. I do not evaluate the three-way N:P:K interaction because blocking reduces degrees of freedom. I determined (and standardized) soil P concentration (P.conc) at five depths. Depth is the repeated measures variable. I retain three depths to illustrate my problem. The data follow:

z = data.frame(plot=c(1:7, 9:15, 17:24, 26:30, 32:36))
z$P = factor(ifelse(z$plot %in% c(1,3,5,7,13:15,17,20:22,24,29,30,33,35), 1, 0))
z$K = factor(ifelse(z$plot %in% c(2,5,7,9,10,13,17:21,27,29,32:34), 1, 0))
z$N = factor(ifelse(z$plot %in% c(3:5,9:11,15,17,20,22,23,27:29,34,35), 1, 0))
z$block = factor(ifelse(z$plot %in% c(1,2,4,5,9,12,13,15,19,20,23,24,28:30,32), 1, 2))
z$rep = ifelse(z$plot %in% c(1:7,10), 1, 0)
z$rep = ifelse(z$plot %in% c(9,11:15,17,18), 2, z$rep)
z$rep = ifelse(z$plot %in% c(19:24,26,27), 3, z$rep)
z$rep = factor(ifelse(z$plot %in% c(28:30,32:36), 4, z$rep))
z$plot = factor(z$plot)
z = rbind(z,z,z)
z$depth = factor(rep(c("0-5 cm", "10-20 cm", "50-100 cm"), each=32))
z$P.conc = c(1.59, -0.34, 1.78, -0.32, 1.76, -0.29, 3.42, -0.44, -0.51, -0.48, -0.56,
             3.55, 1.87, 0.69, 2.40, -0.17, -0.08, 2.03, 1.62, 1.16, -0.59, 1.91, -0.17,
             0.00, -0.46, 1.28, 1.86, -0.12, 1.44, -0.50, 2.21, -0.45, 0.81, -0.72, -0.10,
            -0.53, 0.55, -0.84, 1.95, -0.90, -0.75, -0.81, -0.68, 1.02, 0.63, 0.10, 0.13,
            -0.78, -0.70, 0.34, 0.36, 0.71, -0.86, 1.02, -0.94, -0.26, -0.67, 0.75, 0.48,
            -0.77, 0.51, -0.44, 1.00, -0.81, -0.37, -0.85, -0.72, -0.68, -0.62, -0.60, 0.54,
            -0.94, -0.88, -0.85, -0.77, -0.56, -0.81, -0.66, -0.62, -0.73, -0.73, -0.52, -0.54,
            -0.57, -1.12, -0.97, -0.80, -0.50, -0.98, -0.56, -0.63, -0.77, -0.33, -0.96, -0.83,
            -0.93)

I can use the aov{stats} function to analyze these data as follows:

m.aov = aov(P.conc ~ (N*P + N*K + P*K + rep/block)*depth + Error(plot/depth), data=z)
summary(m.aov)

I can use the lmer{lme4} function to analyze these data as follows:

library(lme4)
m.lmer = lmer(P.conc ~ (N*P + N*K + P*K + rep/block)*depth + (1|plot), data=z)           
anova(m.lmer)

Happily, these analyses with the aov{stats} and lmer{lme4} functions produce identical degrees of freedom, sums of squares, mean squares and F-values. I then used the mauchly.test{stats} function to evaluate the sphericity assumption as follows:

temp = z[1:32, c("N", "P", "K", "rep", "block")]
mdl.mtrx = model.matrix(~ N*P + N*K + P*K + rep/block, temp)
temp = cbind(z[as.character(z$depth)=="0-5 cm","P.conc"],
             z[as.character(z$depth)=="10-20 cm","P.conc"],
             z[as.character(z$depth)=="50-100 cm","P.conc"])
m.lm = lm(temp ~ mdl.mtrx)
mauchly.test(m.lm, X=~1)

Sadly, the sphericity assumption is rejected. The p-value is ~10^-6 when I retain all five depths. So, to access Greenhouse-Geisser corrections, I attempted to implement the anova_test{rstatix} function as follows:

library(rstatix)
anova_test(data=z, formula=P.conc~(N*P+N*K+P*K+rep/block)*depth + Error(plot/depth))

This application of the anova_test{rstatix} function produces the following error message:

Error in solve.default(wcrossprod(model.matrix(mod), w = wts)) : Lapack routine dgesv: system is exactly singular: U[27,27] = 0

So, I next attempted to use the Anova{stats} function as follows:

# m.lm was defined above to implement the mauchly.test{stats} function
ws.factors = unique(z$depth)
m.Anova = Anova(m.lm, idata = data.frame(ws.factors), idesign = ~ws.factors, type=3)

This application of the Anova{stats} function produces the following the following error message:

Error in solve.default(crossprod(model.matrix(mod))) : Lapack routine dgesv: system is exactly singular: U[2,2] = 0

The two error messages are very similar. This is unsurprising because anova_test{rstatix} is a wrapper around Anova{car} and aov{stats}. Both error messages state "... system is exactly singular ...".

But, the analysis works with the aov{stats} function so we know the system is not singular. The problem must lie with my implementation of the anova_test{rstatix} and Anova{car} functions. What have I done wrong?

Answer 1

The coding part of this question is off-topic here, but I suspect that there is a statistical aspect having to do with the different flavors of ANOVA testing.

Your model may be close to overfitting in any event; I count 21 fixed effects to estimate plus several random effects in your model, with only 96 observations. When I see warnings about singularities like those, it has meant that I am trying to ask too much of a limited amount of data.

My guess is that the ability to get a statistical test result with aov() but not with Anova() (the latter with a capital "A" from the car package, not base stats) has to do with the sequential testing done by aov() or anova(), evaluating by default one predictor at a time in the order in which they were entered into the model. The test performed by Anova() works differently and might be more sensitive to such problems.

Also, please do consider the issues with estimating p-values in mixed models. This is not so straightforward as in models without random effects, and there may be better approaches than what you are using even if you solve your singularity problem.