rvegan

How can the effect size of a PERMANOVA be calculated?


I use the "vegan" package to perform a PERMANOVA (adonis2()), and I also want to calculate the effect size (ω²). For this, I tried to use omega_squared() from the "effectsize" package, but I failed. I think it does not understand the output table, specifically the part with the mean squares. Is it possible to fix this or do I have to calculate manually?

library(vegan)
#> Lade nötiges Paket: permute
#> Lade nötiges Paket: lattice
#> This is vegan 2.6-4
library(effectsize)

data(dune)
data(dune.env)

ado <- adonis2(dune ~ Management, data = dune.env, permutations = 100)

ado
#> Permutation test for adonis under reduced model
#> Terms added sequentially (first to last)
#> Permutation: free
#> Number of permutations: 100
#> 
#> adonis2(formula = dune ~ Management, data = dune.env, permutations = 100)
#>            Df SumOfSqs      R2      F   Pr(>F)   
#> Management  3   1.4686 0.34161 2.7672 0.009901 **
#> Residual   16   2.8304 0.65839                   
#> Total      19   4.2990 1.00000                   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

w2 <- omega_squared(ado)
#> Error in `[[<-.data.frame`(`*tmp*`, "Mean_Square", value = numeric(0)): Ersetzung hat 0 Zeilen, Daten haben 3
interpret_omega_squared(w2)
#> Error in interpret(es, rules): Objekt 'w2' nicht gefunden

Created on 2022-11-15 with reprex v2.0.2

EDIT

I tried to do it manually:

library(vegan, quietly = T, warn.conflicts = F)
#> This is vegan 2.6-4
library(effectsize)
library(dplyr, quietly = T, warn.conflicts = F)
library(tibble)
library(purrr)

data(dune)
data(dune.env)

ado <- adonis2(dune ~ Management, data = dune.env, permutations = 100)

w2 <- omega_squared(ado) # Does not work
#> Error in `[[<-.data.frame`(`*tmp*`, "Mean_Square", value = numeric(0)): Ersetzung hat 0 Zeilen, Daten haben 3
interpret_omega_squared(w2) # Does not work
#> Error in interpret(es, rules): Objekt 'w2' nicht gefunden

ado_tidy <- tibble( # manually create Adonis test result table
  parameter = c("Management", "Residual", "Total"),
  df = ado %>% pull("Df"), # Degree of freedom
  ss = ado %>% pull("SumOfSqs"), # sum of squares
  meansqs = ss / df, # mean squares
  p_r2 = ado %>% pull("R2"), # partial R²
  f = ado %>% pull("F"), # F value
  p = ado %>% pull("Pr(>F)") # p value
  )

ado_tidy
#> # A tibble: 3 x 7
#>   parameter     df    ss meansqs  p_r2     f        p
#>   <chr>      <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>
#> 1 Management     3  1.47   0.490 0.342  2.77  0.00990
#> 2 Residual      16  2.83   0.177 0.658 NA    NA      
#> 3 Total         19  4.30   0.226 1     NA    NA

# Formula:
# W2 = (DFm * (F - 1)) / ((DFm * (F - 1)) + (DFm + 1)) 

W2 <- abs(
  (ado_tidy %>% pull(df) %>% chuck(3) * (ado_tidy %>% pull(f) %>% chuck(1) - 1)) /
    ((ado_tidy %>% pull(df) %>% chuck(3) * (ado_tidy %>% pull(f) %>% chuck(1) - 1) +
      ado_tidy %>% pull(df) %>% chuck(3) + 1)
     )
  )

W2
#> [1] 0.6267099
interpret_omega_squared(W2, rules = "field2013")
#> [1] "large"
#> (Rules: field2013)

Created on 2022-11-15 with reprex v2.0.2

Hopefully, the equation is correct...


Solution

  • Here is the MicEco::adonis_OmegaSq function edited so that it works both with the current vegan::adonis2 and deprecated vegan::adonis:

    #' Calculate (partial) Omega-squared (effect-size calculation) for PERMANOVA and add it to the input object
    #'
    #' @param adonisOutput An adonis object
    #' @param partial Should partial omega-squared be calculated (sample size adjusted). Default TRUE
    #' @return Original adonis object with the (partial) Omega-squared values added
    #' @import vegan
    #' @export
    adonis_OmegaSq <- function(adonisOutput, partial = TRUE){
        if(!(is(adonisOutput, "adonis") || is(adonisOutput, "anova.cca")))
            stop("Input should be an adonis object")
        if (is(adonisOutput, "anova.cca")) {
            aov_tab <- adonisOutput
            aov_tab$MeanSqs <- aov_tab$SumOfSqs / aov_tab$Df
            aov_tab$MeanSqs[length(aov_tab$Df)] <- NA
        } else {
            aov_tab <- adonisOutput$aov.tab
        }
        heading <- attr(aov_tab, "heading")
        MS_res <- aov_tab[pmatch("Residual", rownames(aov_tab)), "MeanSqs"]
        SS_tot <- aov_tab[rownames(aov_tab) == "Total", "SumsOfSqs"]
        N <- aov_tab[rownames(aov_tab) == "Total", "Df"] + 1
        if(partial){
            omega <- apply(aov_tab, 1, function(x) (x["Df"]*(x["MeanSqs"]-MS_res))/(x["Df"]*x["MeanSqs"]+(N-x["Df"])*MS_res))
            aov_tab$parOmegaSq <- c(omega[1:(length(omega)-2)], NA, NA)
        } else {
            omega <- apply(aov_tab, 1, function(x) (x["SumsOfSqs"]-x["Df"]*MS_res)/(SS_tot+MS_res))
            aov_tab$OmegaSq <- c(omega[1:(length(omega)-2)], NA, NA)
        }
        if (is(adonisOutput, "adonis"))
            cn_order <- c("Df", "SumsOfSqs", "MeanSqs", "F.Model", "R2",
                          if (partial) "parOmegaSq" else "OmegaSq", "Pr(>F)")
        else
            cn_order <- c("Df", "SumOfSqs", "F", if (partial) "parOmegaSq" else "OmegaSq",
                          "Pr(>F)")
        aov_tab <- aov_tab[, cn_order]
        attr(aov_tab, "names") <- cn_order
        attr(aov_tab, "heading") <- heading
        if (is(adonisOutput, "adonis"))
            adonisOutput$aov.tab <- aov_tab
        else
            adonisOutput <- aov_tab
        return(adonisOutput)
    }
    

    source() this function and it should work. In my test it gave the same results for both adonis2 and adonis.