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

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?

• I get the same error when using Anova() on an lm() object, which has interaction predictor (p1*p2*p3). Making the model with only one predictor at a time all works, but as soon as I add another (either using + or *), the Anova() call will give me > Error in solve.default(crossprod(model.matrix(mod))) : Lapack > routine dgesv: system is exactly singular: U[6,6] = 0 – Teodor Jul 21 at 10:13

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.