16380
2015-05-12 16:35:03 +0000
vcov() fails for maov and manova objects
2015-12-14 13:48:43 +0000
1
1
1
Unclassified
R
Models
R 3.2.0
All
All
CLOSED
FIXED
P5
normal
---
0
russell-lenth
R-core
murdoch
oldest_to_newest
90448
0
russell-lenth
2015-05-12 16:35:03 +0000
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)
}
90449
1
murdoch
2015-05-12 16:46:36 +0000
Could you please expand your example, to make it reproducible? Where did "fit" come from?
90450
2
russell-lenth
2015-05-12 17:09:05 +0000
# 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)
90451
3
murdoch
2015-05-12 17:34:59 +0000
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?
90452
4
russell-lenth
2015-05-12 19:48:59 +0000
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
90453
5
russell-lenth
2015-05-12 19:51:37 +0000
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.
90454
6
murdoch
2015-05-12 21:09:00 +0000
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?
90455
7
russell-lenth
2015-05-12 21:40:34 +0000
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...
90458
8
murdoch
2015-05-13 13:56:51 +0000
I have put the change into R-devel, rev 68361.