Bugzilla – Full Text Bug Listing

Summary: | vcov() fails for maov and manova objects | ||
---|---|---|---|

Product: | R | Reporter: | Russ Lenth <russell-lenth> |

Component: | Models | Assignee: | 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
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. |