Extracting Multivariate Tests from the output of Anova or Manova function from car package

Question

I wonder how to extract the Multivariate Tests: Site portion from the output of fm1 in the following 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

Rich Scriven · Answer 1 · 2014-09-17T23:56:43.953

fm1$multivariate.tests gets you to the Site portion of the fm1 output.

Then you could use a combination of cat and capture.output for nice printing, or just capture.output for a character 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 of

unlist(fm1$multivariate.tests, recursive = FALSE)

it doesn't look like the results are easily accessible as numeric values. So, as you requested, here is what it took to manipulate the results into a matrix. Having done this and then seen user20650's answer, I recommend you follow that suggestion and get the values via an 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 feels like improving my second code chunk, feel free.

(+1): Nice workaround. Can this be transformed to data.frame to use for xtable? Thanks — MYaseen208, Sep 17 '14 at 19:51
Hmm, well I could try to get it into a matrix with some string manipulation — Rich Scriven, Sep 17 '14 at 20:02

user20650 · Accepted Answer · 2014-09-18T00:05:32.733

I also couldn't find how to extract the table of tests but as a workaround you can calculate the results by running the Anova command over all test types.

However the print method, print.Anova.mlm does not return the results, so this needs to be tweaked a little.

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

As the output table is of class anova and data.frame you can use xtable on it.

xtable:::xtable(tab)

+1. This method never even crossed my mind. I do so much data manipulation that I forgot about statistics. — Rich Scriven, Sep 17 '14 at 21:31

Extracting Multivariate Tests from the output of Anova or Manova function from car package

2 Answers2