How to include repeated measures in adonis function (vegan package)?

General
1

Hello everyone,
I am working with forest restoration and I want to test the effect of 3 restoration treatments in the species composition.
To summarize my data:
I have presence/absence and abundance data for several vascular plants species in 18 sites ("ID").
In these sites, 3 different treatments were applied (6 sites within each "Treatment") and the data was collected before, right after and 8 years after ("Timeline") the treatments were applied, therefore I am also working with repeated measures.

Before performing the permanova test with adonis2(), I applied a square root transformation to my data and then created a distance matrix using vegdist(), method bray (FL.dist)

The assumption of homogeneity of dispersion among groups was tested and fulfilled. To perform the permanova I used the following code:
FL.div<-adonis2(FL.dist~Treatment*Timeline, data=Permanova_fieldlayer,
permutations = 999,method="bray", strata="ID")

My question would be, did I treat Timeline as another fixed factor effect? How can I specify that this variable should be use as repeated measures?
Also, I found a new function to peform posthoc with R called pairwise.adonis2 ( by Martinez Arbizu, P. (2019)) has anyone use it ?

Thanks,
Clara Espinosa

2

Hi, and welcome.

A reproducible example, called a reprex would be a huge help in parsing your question. Could you post one?

The examples on the adonis2 help page suggest to me as a non-specialist that FL.div model has one fixed and one interaction term, unless the fixed term is expressly excluded.

library(vegan)
#> Loading required package: permute
#> Loading required package: lattice
#> This is vegan 2.5-6
data(dune)
data(dune.env)
# default test by terms
adonis2(dune ~ Management*A1, data = dune.env)
#> Permutation test for adonis under reduced model
#> Terms added sequentially (first to last)
#> Permutation: free
#> Number of permutations: 999
#> 
#> adonis2(formula = dune ~ Management * A1, data = dune.env)
#>               Df SumOfSqs      R2      F Pr(>F)    
#> Management     3   1.4686 0.34161 3.2629  0.001 ***
#> A1             1   0.4409 0.10256 2.9387  0.019 *  
#> Management:A1  3   0.5892 0.13705 1.3090  0.190    
#> Residual      12   1.8004 0.41878                  
#> Total         19   4.2990 1.00000                  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# overall tests
adonis2(dune ~ Management*A1, data = dune.env, by = NULL)
#> Permutation test for adonis under reduced model
#> Permutation: free
#> Number of permutations: 999
#> 
#> adonis2(formula = dune ~ Management * A1, data = dune.env, by = NULL)
#>          Df SumOfSqs      R2      F Pr(>F)    
#> Model     7   2.4987 0.58122 2.3792  0.001 ***
#> Residual 12   1.8004 0.41878                  
#> Total    19   4.2990 1.00000                  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Created on 2019-11-11 by the reprex package (v0.3.0)

1 Like
3

Hi I tried to create a reprex, can you read it with this format or how should I load it?

reprex%20permanova

1 Like
4

That is not a reprex, that is a screenshot, please follow this guide and try to make a proper reproducible example.

5

Half way there: Just cut and paste the output into your message.

6

Thank you ofr your kind response, I think it is working now

MULTIVARIANCE ANALYSIS - PERMANOVA: Permutational multivariate analysis of variance
library(reprex)
library(vegan)
library(readxl)

load data

Permanova_fieldlayer <- read_excel("C:/Users/Usuario/Dropbox/SLU restoration/Permanova_fieldlayer.xlsx")

str(Permanova_fieldlayer)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 54 obs. of 37 variables:
Treatment : chr "Burning" "Burning" "Burning" "Burning" ... ID : num 1935 1935 1935 2746 2746 ...
Timeline : chr "Before" "After" "Follow up" "Before" ... Vaccinium myrtillus : num 23 24 21 24 23 22 24 21 21 24 ...
Vaccinium vitis-idaea : num 22 24 25 24 21 24 23 22 25 23 ... Deschampsia flexuosa : num 7 7 17 1 2 3 4 5 16 10 ...
Linnaea borealis : num 9 0 2 8 0 3 7 1 3 3 ... Luzula pilosa : num 0 0 0 0 0 0 1 0 0 0 ...
Juniperus comunis : num 0 0 0 0 0 0 0 0 0 0 ... Melampyrum sylvaticum : num 0 0 0 0 0 2 1 0 4 0 ...
Lycopodium annotinum : num 0 0 0 0 0 0 2 0 0 1 ... Maianthemum bifolium : num 0 0 0 0 0 0 0 0 0 0 ...
Empetrum nigrum : num 2 0 0 6 0 0 1 0 0 3 ... Calluna vulgaris : num 0 0 0 1 0 0 0 0 0 0 ...
Andromeda polifolia : num 0 0 0 0 0 0 0 0 0 0 ... Rubus chamaemorus : num 0 0 0 0 0 0 1 0 1 0 ...
Vaccinium uglinosum : num 0 0 0 0 0 0 0 0 0 0 ... Gymnocarpium dryopteris : num 0 0 0 0 0 0 0 0 0 0 ...
Trientalis europaea : num 1 0 0 0 0 0 0 0 0 0 ... Melampyrum pratense : num 0 0 2 1 0 1 0 0 3 0 ...
Epibolium angustifolium : num 0 4 9 0 2 8 0 9 12 0 ... Equisetum pratense : num 0 0 0 0 0 1 0 0 0 0 ...
Equisetum sylvaticum : num 0 0 0 0 1 0 0 0 3 1 ... Orthilia secunda : num 0 0 0 0 0 0 0 0 0 0 ...
Deschampsia cespitosa : num 0 0 0 0 0 0 0 0 1 0 ... Rubus idaeus : num 0 0 1 0 0 0 0 1 0 0 ...
Solidago virgaurea : num 0 0 0 0 0 0 0 0 0 0 ... Calamagrostis purpurea : num 0 0 0 0 0 0 0 0 0 0 ...
Godyera repen : num 0 0 0 0 0 0 0 0 0 0 ... Vaccinium oxycoccos : num 0 0 0 2 0 0 0 0 0 0 ...
Diphasiastrum complanatum: num 0 0 0 0 0 0 0 0 0 0 ... Listera cordata : num 0 0 0 0 0 0 0 0 0 0 ...
Ranunculus lapponicus : num 0 0 0 0 0 0 0 0 0 0 ... Melampyrum sp. : num 0 0 0 0 0 0 0 0 0 0 ...
Carex sp : num 0 0 0 2 0 0 0 0 0 0 ... Rhododendron tomentosum : num 0 0 0 1 0 0 0 0 0 0 ...
$ Equisetum palustre : num 0 0 0 0 0 0 0 0 0 0 ...

Permanova_fieldlayer$Treatment <- as.factor(Permanova_fieldlayer$Treatment)
Permanova_fieldlayer$Timeline <- as.factor(Permanova_fieldlayer$Timeline)
Permanova_fieldlayer$ID <- as.factor(Permanova_fieldlayer$ID)
Permanova_fieldlayer$Group <- paste(Permanova_fieldlayer$Treatment, Permanova_fieldlayer$Timeline)

  1. Transform the data
    FL.mat<-sqrt(FL.matrix)#square root transform

2. Calculate dissimilarity matrix between sites based on species occurence

FL.dist<-vegdist(FL.mat, method='bray')
set.seed(36)

3.test homogeneity of dispersion among groups (different treatments and timeline),

which is a condition (assumption) for adonis.

dispersion<-betadisper(FL.dist, group=Permanova_fieldlayer$Group)

WARNING### some squared distances are negative and changed to zero!!!!!!

permutest(dispersion)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 8 0.028468 0.0035585 0.9517 999 0.479
Residuals 45 0.168254 0.0037390

plot(dispersion, hull=FALSE, ellipse=TRUE) ##sd ellipse

Results show homogeonous dispersion

perform permanova to test if field layer community is influenced by treatment and timeline

FL.div
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = FL.dist ~ Treatment + Timeline, data = Permanova_fieldlayer, permutations = 999, method = "bray", strata = "ID")
Df SumOfSqs R2 F Pr(>F)
Treatment 2 0.27191 0.10486 3.0737 0.001 ***
Timeline 2 0.15391 0.05935 1.7399 0.069 .
Residual 49 2.16737 0.83579
Total 53 2.59319 1.00000

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

does it really take into account repeated measures within each treatment??

POST HOC

from github.com/pmartinezarbizu/pairwiseAdonis####

library(devtools)
install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
install.packages("Rtools")
library(pairwiseAdonis)

Martinez Arbizu, P. (2019). pairwiseAdonis: Pairwise multilevel comparison using adonis.

R package version 0.3

For strata, extract factors into a new data.frame

fac <- data.frame(Treatment=Permanova_fieldlayer$Treatment, Timeline = rep( c('before','after' .... [TRUNCATED]

dim(FL.dist)
NULL

posthoc_FL <-pairwise.adonis2(FL.dist~Treatment+Timeline, data=fac, strata="ID")
Error in [.data.frame(mdat1, , strata) : undefined columns selected

In addition: Warning messages:
1: In betadisper(FL.dist, group = Permanova_fieldlayer$Group) :
some squared distances are negative and changed to zero
2: package ‘devtools’ was built under R version 3.5.3
3: package ‘usethis’ was built under R version 3.5.3

source('~/.active-rstudio-document', echo=TRUE)

7

Great. Please mark your answer as the solution (no false modest!). It will help any others with a similar question.

8

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.