R: Need example of performing boxM() with two-way (factorial) MANOVA - I'm getting an error

Question

I am trying to run Box's M-test for Homogeneity of Covariance Matrices, for two-way MANOVA.

I have conducted searches for an example since yesterday afternoon. I see many examples of using boxM with one-way MANOVA. In every case, if the source also covers two-way MANOVA, they don't include demonstrating running the boxM test in the two-way case. I just need a working example. Once I have the syntax, I'll be able to get this to work.

The boxM function in the biotools package says it is for one classification factor (one-way MANOVA).
https://www.rdocumentation.org/packages/biotools/versions/3.1/topics/boxM

The boxM function in the heplots package says it works with one or more classification factors --
https://www.rdocumentation.org/packages/heplots/versions/1.3-5/topics/boxM

-- However, I get an error when I try to use it: "Model must be completely crossed formula only."

Below, I show that I don't get an error when using either factor alone, but any arrangement of crossing the factors give this error. Note: I don't get this error running Levene's test with the variables crossed.

Response1, Response2, and Response3 are continuous.
Factor1 has 2 levels. Factor2 has 5 levels.

library(heplots)

> Model2 <- lm(cbind(Response1, Response2, Response3) ~ Factor1, data=Data40)
> boxM(Model2)

    Box's M-test for Homogeneity of Covariance Matrices

data:  Y
Chi-Sq (approx.) = 3.5562, df = 6, p-value = 0.7365

> Model2 <- lm(cbind(Response1, Response2, Response3) ~ Factor2, data=Data40)
> boxM(Model2)

    Box's M-test for Homogeneity of Covariance Matrices

data:  Y
Chi-Sq (approx.) = 35.079, df = 24, p-value = 0.06724

> Model2 <- lm(cbind(Response1, Response2, Response3) ~ Factor1 * Factor2, data=Data40)
> boxM(Model2)
Error in boxM.formula(formula(Y), data = eval(data, envir = environment(formula(Y))),  : 
  Model must be completely crossed formula only.

> Model2 <- lm(cbind(Response1, Response2, Response3) ~ Factor1 + Factor2 + Factor1 * Factor2, data=Data40)
> boxM(Model2)
Error in boxM.formula(formula(Y), data = eval(data, envir = environment(formula(Y))),  : 
  Model must be completely crossed formula only.

> Model2 <- lm(cbind(Response1, Response2, Response3) ~ Factor1 + Factor2 + Factor1:Factor2, data=Data40)
> boxM(Model2)
Error in boxM.formula(formula(Y), data = eval(data, envir = environment(formula(Y))),  : 
  Model must be completely crossed formula only.

Answer 1

Don't know the package never used it but in a couple minutes of sleuthing it appears you may be specifying the formula in a way it doesn't like... Using iris since the package author does and you provide no data.

library(heplots)

# adding a bogus second factor to iris
iris$nonsense <- rep(1:2)
iris$nonsense <- factor(iris$nonsense)

# one factor
boxM( cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ nonsense, data=iris)
#> 
#>  Box's M-test for Homogeneity of Covariance Matrices
#> 
#> data:  Y
#> Chi-Sq (approx.) = 16.389, df = 10, p-value = 0.08904

# second factor
boxM( cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, data=iris)
#> 
#>  Box's M-test for Homogeneity of Covariance Matrices
#> 
#> data:  Y
#> Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16

# crossed note not including the `lm`

boxM( cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species * nonsense, data=iris)
#> 
#>  Box's M-test for Homogeneity of Covariance Matrices
#> 
#> data:  Y
#> Chi-Sq (approx.) = 169.1, df = 50, p-value = 7.609e-15