Bug 16380

Summary: vcov() fails for maov and manova objects
Product: R Reporter: Russ Lenth <russell-lenth>
Component: ModelsAssignee: R-core <R-core>
Status: CLOSED FIXED    
Severity: normal CC: murdoch
Priority: P5    
Version: R 3.2.0   
Hardware: All   
OS: All   

Description Russ Lenth 2015-05-12 16:35:03 UTC
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)
    }
Comment 1 Duncan Murdoch 2015-05-12 16:46:36 UTC
Could you please expand your example, to make it reproducible?  Where did "fit" come from?
Comment 2 Russ Lenth 2015-05-12 17:09:05 UTC
# 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)
Comment 3 Duncan Murdoch 2015-05-12 17:34:59 UTC
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?
Comment 4 Russ Lenth 2015-05-12 19:48:59 UTC
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
Comment 5 Russ Lenth 2015-05-12 19:51:37 UTC
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.
Comment 6 Duncan Murdoch 2015-05-12 21:09:00 UTC
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?
Comment 7 Russ Lenth 2015-05-12 21:40:34 UTC
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...
Comment 8 Duncan Murdoch 2015-05-13 13:56:51 UTC
I have put the change into R-devel, rev 68361.