Extract multivariate tests from Anova or Manova function results from a car package
Wondering how to extract Multivariate Tests: Site
some Multivariate Tests: Site
of the fm1 output in the next MWE.
library(car)
fm1 <- summary(Anova(lm(cbind(Al, Fe, Mg, Ca, Na) ~ Site, data=Pottery)))
fm1
Type II MANOVA Tests:
Sum of squares and products for error:
Al Fe Mg Ca Na
Al 48.2881429 7.08007143 0.60801429 0.10647143 0.58895714
Fe 7.0800714 10.95084571 0.52705714 -0.15519429 0.06675857
Mg 0.6080143 0.52705714 15.42961143 0.43537714 0.02761571
Ca 0.1064714 -0.15519429 0.43537714 0.05148571 0.01007857
Na 0.5889571 0.06675857 0.02761571 0.01007857 0.19929286
------------------------------------------
Term: Site
Sum of squares and products for the hypothesis:
Al Fe Mg Ca Na
Al 175.610319 -149.295533 -130.809707 -5.8891637 -5.3722648
Fe -149.295533 134.221616 117.745035 4.8217866 5.3259491
Mg -130.809707 117.745035 103.350527 4.2091613 4.7105458
Ca -5.889164 4.821787 4.209161 0.2047027 0.1547830
Na -5.372265 5.325949 4.710546 0.1547830 0.2582456
Multivariate Tests: Site
Df test stat approx F num Df den Df Pr(>F)
Pillai 3 1.55394 4.29839 15 60.00000 2.4129e-05 ***
Wilks 3 0.01230 13.08854 15 50.09147 1.8404e-12 ***
Hotelling-Lawley 3 35.43875 39.37639 15 50.00000 < 2.22e-16 ***
Roy 3 34.16111 136.64446 5 20.00000 9.4435e-15 ***
---
Signif. codes: 0 ‘*** 0.001 ‘** 0.01 ‘* 0.05 ‘. 0.1 ‘ 1
source to share
I also couldn't find how to retrieve the test table, but as a workaround, you can compute the results by running the command Anova
on all test types.
However, the print method print.Anova.mlm
returns no results, so this needs to be tweaked a bit.
library(car)
# create new print function
outtests <- car:::print.Anova.mlm
# allow the function to return the results and disable print
body(outtests)[[16]] <- quote(invisible(tests))
body(outtests)[[15]] <- NULL
# Now run the regression
mod <- lm(cbind(Al, Fe, Mg, Ca, Na) ~ Site, data=Pottery)
# Run the Anova over all tests
tab <- lapply(c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"),
function(i) outtests(Anova(mod, test.statistic=i)))
tab <- do.call(rbind, tab)
row.names(tab) <- c("Pillai", "Wilks", "Hotelling-Lawley", "Roy")
tab
# Type II MANOVA Tests: Pillai test statistic
# Df test stat approx F num Df den Df Pr(>F)
#Pillai 3 1.554 4.298 15 60.000 2.413e-05 ***
#Wilks 3 0.012 13.089 15 50.091 1.840e-12 ***
#Hotelling-Lawley 3 35.439 39.376 15 50.000 < 2.2e-16 ***
#Roy 3 34.161 136.644 5 20.000 9.444e-15 ***
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Since the output table has a class Anova
and data.frame
, you can use xtable
.
xtable:::xtable(tab)
source to share
fm1$multivariate.tests
takes you to the Site
output part fm1
.
Then you can use a combination cat
and capture.output
for a beautiful print or just capture.output
for the symbol vector.
> cat(capture.output(fm1$multivariate.tests)[18:26], sep = "\n")
#
# Multivariate Tests: Site
# Df test stat approx F num Df den Df Pr(>F)
# Pillai 3 1.55394 4.29839 15 60.00000 2.4129e-05 ***
# Wilks 3 0.01230 13.08854 15 50.09147 1.8404e-12 ***
# Hotelling-Lawley 3 35.43875 39.37639 15 50.00000 < 2.22e-16 ***
# Roy 3 34.16111 136.64446 5 20.00000 9.4435e-15 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Update: From the result
unlist(fm1$multivariate.tests, recursive = FALSE)
it doesn't look like the results are readily available as numeric values. So, as you asked, here's what it takes to manipulate the results into a matrix. After doing this and then seeing user20650's answer, I recommend that you follow this suggestion and get the values through the ANOVA table.
co <- capture.output(fm1$multivariate.tests)[20:24]
s <- strsplit(gsub("([*]+$)|[<]", "", co[-1]), "\\s+")
dc <- do.call(rbind, lapply(s, function(x) as.numeric(x[-1])))
row.names(dc) <- sapply(s, "[", 1)
s2 <- strsplit(co[1], " ")[[1]]
s2 <- s2[nzchar(s2)]
s3 <- s2[-c(1, length(s2))]
colnames(dc) <- c(s2[1], paste(s3[c(TRUE, FALSE)], s3[c(FALSE, TRUE)]), s2[10])
dc
# Df test stat approx F num Df den Df Pr(>F)
# Pillai 3 1.55394 4.29839 15 60.00000 2.4129e-05
# Wilks 3 0.01230 13.08854 15 50.09147 1.8404e-12
# Hotelling-Lawley 3 35.43875 39.37639 15 50.00000 2.2200e-16
# Roy 3 34.16111 136.64446 5 20.00000 9.4435e-15
If anyone wants to improve my second piece of code, don't hesitate.
source to share