Do groups differ in composition?
Does community structure vary among regions or over time?
Do environmental variables explain community patterns? Which species are responsible for differences among groups?
## Running Code
library(ggvegan)Loading required package: veganWarning: package 'vegan' was built under R version 4.1.3Loading required package: permuteWarning: package 'permute' was built under R version 4.1.3Loading required package: latticeThis is vegan 2.6-2Loading required package: ggplot2library(vegan) You can add options to executable code like this
birds<-read.csv('https://raw.githubusercontent.com/collnell/lab-demo/master/bird_by_fg.csv')
trees<-read.csv('https://raw.githubusercontent.com/collnell/lab-demo/master/tree_comp.csv')Birds abundance data: 32 rows means plots 12 of which have 1 tree species with diversity M and 20 with 4 tree species diversity P columns are
head(birds) DIVERSITY PLOT CA FR GR HE IN NE OM
1 M 3 0 0 0 0 2 0 0
2 M 9 0 0 2 0 6 0 4
3 M 12 0 0 0 0 2 0 2
4 M 17 0 0 0 0 7 0 4
5 M 20 0 0 0 0 1 0 4
6 M 21 0 0 3 0 14 0 7unique(birds$DIVERSITY)[1] "M" "P"Bird abundances are totalled according to their feeding guild (columns).
colnames(birds)[1] "DIVERSITY" "PLOT" "CA" "FR" "GR" "HE"
[7] "IN" "NE" "OM" columns are species 6 species A to F
head(trees) PLOT comp A B C D E F row col
1 3 D 0 0 0 1 0 0 3 1
2 9 A 1 0 0 0 0 0 2 2
3 12 E 0 0 0 0 1 0 5 2
4 17 F 0 0 0 0 0 1 3 3
5 20 A 1 0 0 0 0 0 6 3
6 21 B 0 1 0 0 0 0 7 3Question: Is C. pentandara (B) associated with variation in bird species composition? Or D & F (both Fabaceae)? Parametric test for differences between independent groups for multiple continuous dependent variables. Like ANOVA for many response variable
bird.matrix<- as.matrix(birds[,3:9]) # response variablestrees$B<-as.factor(trees$B)birds.manova<-manova(bird.matrix~as.factor(B), data=trees)
summary(birds.manova) Df Pillai approx F num Df den Df Pr(>F)
as.factor(B) 1 0.39147 2.2056 7 24 0.07027 .
Residuals 30
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary.aov(birds.manova) Response CA :
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(B) 1 0.2449 0.24494 1.3983 0.2463
Residuals 30 5.2551 0.17517
Response FR :
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(B) 1 0.0062 0.006199 0.0686 0.7952
Residuals 30 2.7125 0.090418
Response GR :
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(B) 1 5.686 5.6862 2.6732 0.1125
Residuals 30 63.814 2.1271
Response HE :
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(B) 1 0.02138 0.021382 0.6771 0.4171
Residuals 30 0.94737 0.031579
Response IN :
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(B) 1 73.65 73.655 2.7977 0.1048
Residuals 30 789.81 26.327
Response NE :
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(B) 1 0.3968 0.39676 1.5655 0.2205
Residuals 30 7.6032 0.25344
Response OM :
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(B) 1 19.44 19.441 1.4069 0.2449
Residuals 30 414.56 13.819 Problem: Most ecological data is overdispersed, has many 0’s or rare species, unequal sample sizes. Solution: Dissimilarity coefficients, permutation tests
PERMANOVA: Permutational multivariate analysis of variance Non-parametric, based on dissimilarities. Allows for partitioning of variability, similar to ANOVA, allowing for complex design (multiple factors, nested design, interactions, covariates).
Uses permutation to compute F-statistic (pseudo-F). Null hypothesis: Groups do not differ in spread or position multivariate space.
PERMANOVA (which is basically adonis()) was found to be largely unaffected by heterogeneity in Anderson & Walsh’s simulations but only for balanced designs.
For unbalanced designs PERMANOVA and ANOSIM were
1.too liberal if the smaller group had greater dispersion, and
2. too conservative if the larger group had greater dispersion.This result was especially so for ANOSIM.
# sqrt transformation
bird.mat<- sqrt(bird.matrix)Quantify pairwise compositional dissimilarity between sites based on species occurrences.
Bray-Curtis dissimilarity (abundance weighted)
Jaccard (presence/absence)
Gower’s (non-continuous variables)
Dissimilarity: 0 = sites are identical, 1 = sites do not share any species ## When to use :
bird.dist<-vegdist(bird.mat, method='bray')
bird.dist 1 2 3 4 5 6 7
2 0.61136859
3 0.33333333 0.34919819
4 0.53326015 0.15324060 0.24314702
5 0.54691816 0.32308204 0.17157288 0.21525044
6 0.70332241 0.16131945 0.48329258 0.27212306 0.46040543
7 0.33333333 0.34919819 0.00000000 0.24314702 0.17157288 0.48329258
8 0.56548939 0.18797110 0.28608138 0.04614414 0.25882380 0.22885259 0.28608138
9 0.31783725 0.36432551 0.13165250 0.25938626 0.30216948 0.49646626 0.13165250
10 0.64031869 0.27728093 0.39028832 0.16256866 0.36504508 0.15009010 0.39028832
11 0.00000000 0.61136859 0.33333333 0.53326015 0.54691816 0.70332241 0.33333333
12 0.54691816 0.32308204 0.17157288 0.21525044 0.00000000 0.46040543 0.17157288
13 0.61401441 0.34387252 0.67645606 0.49887954 0.77562355 0.25944081 0.67645606
14 0.58356587 0.25605328 0.31066339 0.12649698 0.34778817 0.35130514 0.31066339
15 0.56548939 0.18797110 0.28608138 0.04614414 0.25882380 0.22885259 0.28608138
16 0.72754123 0.23802437 0.52041591 0.37762807 0.49861728 0.26303429 0.52041591
17 0.61324320 0.24345665 0.35182873 0.15607334 0.32576539 0.30106633 0.35182873
18 0.54691816 0.32308204 0.17157288 0.21525044 0.00000000 0.46040543 0.17157288
19 0.71599127 0.38237706 0.50261358 0.29559037 0.48028238 0.23339705 0.50261358
20 0.60837273 0.23756138 0.34500284 0.11118252 0.31880333 0.16596176 0.34500284
21 0.65143042 0.14320490 0.40632872 0.30052063 0.43670040 0.23072709 0.40632872
22 0.66463904 0.08973217 0.42559464 0.20350694 0.40118355 0.07263877 0.42559464
23 0.71681417 0.18770497 0.50387598 0.29713243 0.48158190 0.09260034 0.50387598
24 0.68796639 0.22508897 0.46015731 0.31098623 0.43663461 0.33593482 0.46015731
25 0.55683487 0.12626874 0.27444097 0.22361114 0.24699918 0.24078868 0.27444097
26 0.65485688 0.29621346 0.41130566 0.18684455 0.38654863 0.17989575 0.41130566
27 0.41421356 0.74566690 0.58578644 0.69905585 0.60000000 0.80236099 0.58578644
28 0.64731795 0.28632627 0.40037448 0.17418427 0.37536120 0.23317526 0.40037448
29 0.72136194 0.25351449 0.51086912 0.33502554 0.48878247 0.27679305 0.51086912
30 0.22514823 0.44786886 0.44151844 0.35015208 0.61803399 0.56814021 0.44151844
31 0.66761358 0.31329441 0.42996338 0.20862715 0.40566055 0.33022132 0.42996338
32 0.62589308 0.24055117 0.36968817 0.41692739 0.47208752 0.31607884 0.36968817
8 9 10 11 12 13 14
2
3
4
5
6
7
8
9 0.30191672
10 0.11730448 0.40488183
11 0.56548939 0.31783725 0.64031869
12 0.25882380 0.30216948 0.36504508 0.54691816
13 0.51934025 0.59932397 0.43960545 0.61401441 0.77562355
14 0.16398662 0.32623912 0.25971183 0.58356587 0.34778817 0.53136850
15 0.00000000 0.30191672 0.11730448 0.56548939 0.25882380 0.51934025 0.16398662
16 0.33359286 0.53294114 0.39213319 0.72754123 0.49861728 0.50938717 0.30732384
17 0.10880822 0.36692360 0.20654465 0.61324320 0.32576539 0.58527407 0.08102494
18 0.25882380 0.30216948 0.36504508 0.54691816 0.00000000 0.77562355 0.34778817
19 0.25289566 0.51545558 0.13973658 0.71599127 0.48028238 0.52081554 0.22746785
20 0.06537377 0.36018143 0.08435507 0.60837273 0.31880333 0.46043677 0.21725939
21 0.32717581 0.42069827 0.39646820 0.65143042 0.43670040 0.37082702 0.34291636
22 0.23305835 0.43968366 0.16765234 0.66463904 0.40118355 0.24562270 0.29376431
23 0.28570473 0.51669587 0.22231816 0.71681417 0.48158190 0.28265912 0.22907030
24 0.33528532 0.47371155 0.39912955 0.68796639 0.43663461 0.52284012 0.23732163
25 0.21137297 0.29039204 0.30491245 0.55683487 0.24699918 0.49786370 0.33034446
26 0.14192405 0.42560395 0.10604766 0.65485688 0.38654863 0.45422800 0.27987962
27 0.71812092 0.57735027 0.76329932 0.41421356 0.60000000 0.74727025 0.45783311
28 0.12907760 0.41482818 0.16515165 0.64731795 0.37536120 0.52073395 0.26935275
29 0.29164183 0.52356555 0.35481067 0.72136194 0.48878247 0.52908883 0.40748254
30 0.38999489 0.30273374 0.48510664 0.22514823 0.61803399 0.45124872 0.41263301
31 0.16406242 0.44398704 0.24764233 0.66761358 0.40566055 0.59323770 0.29803785
32 0.44024145 0.38455647 0.50042461 0.62589308 0.47208752 0.38107870 0.45396629
15 16 17 18 19 20 21
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 0.33359286
17 0.10880822 0.23508659
18 0.25882380 0.49861728 0.32576539
19 0.25289566 0.32823879 0.18317671 0.48028238
20 0.06537377 0.33917554 0.16306894 0.31880333 0.19067426
21 0.32717581 0.18525545 0.37009734 0.43670040 0.47941690 0.36551714
22 0.23305835 0.29312594 0.31115755 0.40118355 0.27970023 0.19514178 0.21539970
23 0.28570473 0.21480921 0.21430313 0.48158190 0.20067075 0.24927568 0.25360358
24 0.33528532 0.24944365 0.22710649 0.43663461 0.35337277 0.37049038 0.30760103
25 0.21137297 0.28695948 0.26970040 0.24699918 0.41268956 0.26353668 0.23483180
26 0.14192405 0.29479016 0.22727238 0.38654863 0.20265310 0.07726716 0.41129747
27 0.71812092 0.63526633 0.49360547 0.60000000 0.62066362 0.74385341 0.77011942
28 0.12907760 0.28685723 0.21644729 0.37536120 0.25470017 0.09039114 0.40354697
29 0.29164183 0.25973301 0.33054233 0.48878247 0.41032140 0.29936876 0.35633024
30 0.38999489 0.60080816 0.45026308 0.61803399 0.58517583 0.44404748 0.49955046
31 0.16406242 0.29205509 0.24597563 0.40566055 0.34843348 0.21024034 0.27974325
32 0.44024145 0.39648959 0.47759225 0.47208752 0.57167444 0.47361767 0.24051839
22 23 24 25 26 27 28
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 0.09965125
24 0.28615523 0.23125165
25 0.21052774 0.26625369 0.24997612
26 0.20812940 0.25751955 0.41292140 0.32450914
27 0.77826473 0.62170183 0.58550432 0.71297849 0.77222841
28 0.26607207 0.31010058 0.40570740 0.31428514 0.05625373 0.76759188
29 0.30828464 0.34310334 0.36094732 0.27498149 0.25682013 0.81372279 0.21067043
30 0.51683052 0.58628637 0.54764537 0.51915438 0.50402152 0.52786405 0.49419492
31 0.34172875 0.37687830 0.42544349 0.34212046 0.26560205 0.78010478 0.25621221
32 0.30723563 0.33718249 0.39894168 0.25389213 0.51322556 0.75448869 0.50653860
29 30 31
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30 0.59243270
31 0.25479959 0.52073863
32 0.30843495 0.58673774 0.37374728perMANOVA Do monoculture and polyculture plots differ in feeding guild composition?
set.seed(36) #reproducible results
bird.div<-adonis2(bird.dist~as.factor(birds$DIVERSITY), data=birds, permutations=9999)
bird.divPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 9999
adonis2(formula = bird.dist ~ as.factor(birds$DIVERSITY), data = birds, permutations = 9999)
Df SumOfSqs R2 F Pr(>F)
as.factor(birds$DIVERSITY) 1 0.32857 0.12174 4.1585 0.0077 **
Residual 30 2.37033 0.87826
Total 31 2.69890 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Betadisper tests whether two or more groups (for example, grazed and ungrazed sites) are homogeneously dispersed in relation to their species in studied samples. This test can be done to see if one group has more compositional variance than another. betadisper() is an implementation of Marti Anderson’s Permdisp method.
dispersion<-betadisper(bird.dist, group=birds$DIVERSITY)
dispersion
Homogeneity of multivariate dispersions
Call: betadisper(d = bird.dist, group = birds$DIVERSITY)
No. of Positive Eigenvalues: 13
No. of Negative Eigenvalues: 13
Average distance to median:
M P
0.2309 0.2531
Eigenvalues for PCoA axes:
(Showing 8 of 26 eigenvalues)
PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
1.28713 0.69577 0.46833 0.37989 0.17181 0.11120 0.07211 0.05095 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 1 0.00369 0.0036924 0.2231 999 0.64
Residuals 30 0.49659 0.0165530 anova(dispersion)Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 1 0.00369 0.0036924 0.2231 0.6401
Residuals 30 0.49659 0.0165530 plot(dispersion, hull=FALSE, ellipse=TRUE) ##sd ellipse
birdMDS<-metaMDS(bird.mat, distance="bray", k=2, trymax=35, autotransform=TRUE) ##k is the number of dimensionsRun 0 stress 0.1569259
Run 1 stress 0.1572399
... Procrustes: rmse 0.01130213 max resid 0.05306814
Run 2 stress 0.1572399
... Procrustes: rmse 0.01130502 max resid 0.05307926
Run 3 stress 0.1572399
... Procrustes: rmse 0.0113213 max resid 0.05317232
Run 4 stress 0.3886078
Run 5 stress 0.1569259
... Procrustes: rmse 7.862226e-06 max resid 2.095807e-05
... Similar to previous best
Run 6 stress 0.18157
Run 7 stress 0.170243
Run 8 stress 0.158464
Run 9 stress 0.1671135
Run 10 stress 0.1584574
Run 11 stress 0.1703794
Run 12 stress 0.1627566
Run 13 stress 0.1785574
Run 14 stress 0.171026
Run 15 stress 0.1685473
Run 16 stress 0.1673167
Run 17 stress 0.1572399
... Procrustes: rmse 0.01130537 max resid 0.05308626
Run 18 stress 0.1572399
... Procrustes: rmse 0.01130017 max resid 0.05305432
Run 19 stress 0.1671136
Run 20 stress 0.1569259
... New best solution
... Procrustes: rmse 4.837451e-06 max resid 1.545126e-05
... Similar to previous best
*** Solution reachedbirdMDS ##metaMDS takes eaither a distance matrix or your community matrix (then requires method for 'distance=')
Call:
metaMDS(comm = bird.mat, distance = "bray", k = 2, trymax = 35, autotransform = TRUE)
global Multidimensional Scaling using monoMDS
Data: bird.mat
Distance: bray
Dimensions: 2
Stress: 0.1569259
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'bird.mat' stressplot(birdMDS)
library(ggplot2)NMDS1<-birdMDS$points[,1]
NMDS2<-birdMDS$points[,2]
bird.plot<- cbind(birds,NMDS1, NMDS2,trees)head(bird.plot) DIVERSITY PLOT CA FR GR HE IN NE OM NMDS1 NMDS2 PLOT comp A B C D
1 M 3 0 0 0 0 2 0 0 -0.92785776 0.01829522 3 D 0 0 0 1
2 M 9 0 0 2 0 6 0 4 0.14587417 -0.13372296 9 A 1 0 0 0
3 M 12 0 0 0 0 2 0 2 -0.35737841 0.19813157 12 E 0 0 0 0
4 M 17 0 0 0 0 7 0 4 0.01217651 0.12947460 17 F 0 0 0 0
5 M 20 0 0 0 0 1 0 4 -0.18783764 0.42434011 20 A 1 0 0 0
6 M 21 0 0 3 0 14 0 7 0.37943972 -0.13919083 21 B 0 1 0 0
E F row col
1 0 0 3 1
2 0 0 2 2
3 1 0 5 2
4 0 1 3 3
5 0 0 6 3
6 0 0 7 3# remove duplicated columns
bird.plot <- bird.plot[,!duplicated(colnames(bird.plot))]plot ordination
p<-ggplot(bird.plot, aes(NMDS1, NMDS2, color=DIVERSITY))+
geom_point(position=position_jitter(.1), shape=3)+##separates overlapping points
stat_ellipse(type='t',size =1)+ ##draws 95% confidence interval ellipses
theme_minimal()
p
#plot ordination
p<-ggplot(bird.plot, aes(NMDS1, NMDS2, color=DIVERSITY))+
stat_ellipse(type='t',size =1)+
theme_minimal()+geom_text(data=bird.plot,aes(NMDS1, NMDS2, label=comp), position=position_jitter(.35))+
annotate("text", x=min(NMDS1), y=min(NMDS2), label=paste('Stress =',round(birdMDS$stress,3))) #add stress to plot
p
Fit vector to ordination
fit<-envfit(birdMDS, bird.mat)fit
***VECTORS
NMDS1 NMDS2 r2 Pr(>r)
CA 0.94208 -0.33540 0.0791 0.276
FR 0.94830 -0.31736 0.0780 0.321
GR 0.39403 -0.91910 0.7525 0.001 ***
HE 0.13945 -0.99023 0.0487 0.455
IN 0.97680 -0.21416 0.5083 0.001 ***
NE 0.25813 -0.96611 0.0627 0.408
OM 0.78680 0.61721 0.8193 0.001 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Permutation: free
Number of permutations: 999autoplot(fit)
arrows<- data.frame(fit$vector$arrows,R=fit$vectors$r,P=fit$vectors$pvals)arrows NMDS1 NMDS2 R P
CA 0.9420757 -0.3354004 0.07914145 0.276
FR 0.9483046 -0.3173615 0.07799746 0.321
GR 0.3940310 -0.9190972 0.75254539 0.001
HE 0.1394546 -0.9902285 0.04869158 0.455
IN 0.9767994 -0.2141562 0.50831698 0.001
NE 0.2581308 -0.9661100 0.06273997 0.408
OM 0.7868022 0.6172053 0.81929614 0.001arrows$FG <- rownames(arrows) # filter P less than 0.05
arrows.p<-arrows[arrows$P<0.05,]
arrows.p NMDS1 NMDS2 R P FG
GR 0.3940310 -0.9190972 0.7525454 0.001 GR
IN 0.9767994 -0.2141562 0.5083170 0.001 IN
OM 0.7868022 0.6172053 0.8192961 0.001 OMarrows.p NMDS1 NMDS2 R P FG
GR 0.3940310 -0.9190972 0.7525454 0.001 GR
IN 0.9767994 -0.2141562 0.5083170 0.001 IN
OM 0.7868022 0.6172053 0.8192961 0.001 OMp<-ggplot(data=bird.plot, aes(NMDS1, NMDS2))+
geom_point(data=bird.plot, aes(NMDS1, NMDS2, color=DIVERSITY),position=position_jitter(.1))+##separates overlapping points
stat_ellipse(aes(fill=DIVERSITY), alpha=.2,type='t',size =1, geom="polygon")+ ##changes shading on ellipses
theme_minimal()+
geom_segment(data=arrows.p, aes(x=0, y=0, xend=NMDS1, yend=NMDS2, label=FG, lty=FG), arrow=arrow(length=unit(.2, "cm")*arrows.p$R)) ##add arrows (scaled by R-squared value)Warning: Ignoring unknown aesthetics: labelp
ggplot(data=bird.plot, aes(NMDS1, NMDS2))+
geom_point(data=bird.plot, aes(NMDS1, NMDS2, color=DIVERSITY),position=position_jitter(.1))+##separates overlapping points
stat_ellipse(aes(fill=DIVERSITY), alpha=.2,type='t',size =1, geom="polygon")+ ##changes shading on ellipses
theme_minimal()+
geom_segment(data=arrows, aes(x=0, y=0, xend=NMDS1, yend=NMDS2, label=FG, lty=FG), arrow=arrow(length=unit(.2, "cm")*arrows$R)) ##add arrows (scaled by R-squared value)Warning: Ignoring unknown aesthetics: label
Fixed and random Cross and nested permanova
Principal Coordinates Analysis
CO the interest is primarily in the similarity among the objects and the individual data variables are of secondary importance. For this reason, the aim of PCO is often described as reduction in dimensionality, while retaining the important information about relationships among a set of objects.
BACI DISTLM CAP