The `vcov` method produces an error when used with a `maov` or `manova` object: > class(fit) [1] "manova" "maov" "aov" "mlm" "lm" > vcov(fit) Error in seq_len(ny) : argument must be coercible to non-negative integer In addition: Warning message: In seq_len(ny) : first element used of 'length.out' argument Note that the `vcov.mlm` is used, due to inheritance, and it in turn calls `summary.mlm`. But the culprit is actually `coef.aov`: > stats:::coef.aov function (object, ...) { z <- object$coefficients z[!is.na(z)] } <bytecode: 0x00000000123b44c8> <environment: namespace:stats> When the coefficients comprise a matrix, the removal of missings changes it to a vector, and this messes-up `summary.mlm`. I have found that the following simple method fixes the problem: vcov.maov <- function(object, ...) { class(object) <- c("mlm", "lm") vcov(object) }
Could you please expand your example, to make it reproducible? Where did "fit" come from?
# Here is the code that produced `fit` in the original report: tear <- c(6.5, 6.2, 5.8, 6.5, 6.5, 6.9, 7.2, 6.9, 6.1, 6.3, 6.7, 6.6, 7.2, 7.1, 6.8, 7.1, 7.0, 7.2, 7.5, 7.6) gloss <- c(9.5, 9.9, 9.6, 9.6, 9.2, 9.1, 10.0, 9.9, 9.5, 9.4, 9.1, 9.3, 8.3, 8.4, 8.5, 9.2, 8.8, 9.7, 10.1, 9.2) opacity <- c(4.4, 6.4, 3.0, 4.1, 0.8, 5.7, 2.0, 3.9, 1.9, 5.7, 2.8, 4.1, 3.8,1.6, 3.4, 8.4, 5.2, 6.9, 2.7, 1.9) Y <- cbind(tear, gloss, opacity) rate <- factor(gl(2,10), labels = c("Low", "High")) fit <- manova(Y ~ rate)
Thanks. I don't really like manipulating the class the way you did: it works here, but what if the original class vector was different? Something that appears to work here, but probably has bad effects elsewhere, is just to define a coef.maov method that doesn't mess up the matrix structure of fit$coefficients. For your example this works: coef.maov <- function(object, ...) object$coefficients However, if NA coefficients are a possibility here, we'd need to do something smarter. Do you (or someone else) know if manova() can return NA coefficients?
Yes, it is possible: > set.seed(123) > x1 <- rnorm(length(tear)) > x2 <- 4 - 3*x1 > fit2 <- manova(Y ~ x1 + x2 + rate) > fit2$coefficients tear gloss opacity (Intercept) 6.4932183 9.57560322 3.78396921 x1 -0.0431266 -0.07508433 0.08081397 x2 NA NA NA rateHigh 0.5957788 -0.49993898 0.27917123 > coef(fit2) [1] 6.49321835 -0.04312660 0.59577880 9.57560322 -0.07508433 [6] -0.49993898 3.78396921 0.08081397 0.27917123
Oh, I mean to also say that `NA`s are the result of collinearity in the model matrix, so it should be that `NA`s can occur only in whole rows of the coefficient matrix.
Thanks for the example. The obvious variation on my change (deleting the rows of coefficients that include NA) appears not to work: it gets the names wrong. But just using the original, i.e. coef.maov <- function(object, ...) object$coefficients seems okay. Do you see any other problems with it?
I don't see anyproblems with it, based on making that change and testing `coef()`, `vcov()`, `summary()`, and `summary.aov()` on both `fit` and `fit2` (above examples). I also tested `anova(fit, fit2)`, `example(manova)`, and `example(summary.manova)`. That said, I can easily imagine there could be something else in **stats**, or something in one or more contributed packages, that might rely on `coef()` returning a vector, or assumes all results are non-missing. At least I am fairly confident in saying that things work as documented...
I have put the change into R-devel, rev 68361.