rstatisticsanovamultivariate-testingmanova

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


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

Solution

  • 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)