Assignment 14 PCA
Load libraries and data
library("here")
library("vegan")
library("tidyverse")
zp <- read.csv(here("Week_14", "Data", "mangrove_fishpond_zooplankton_edit.csv"))PCA
- Start by characterizing how the environmental variables (columns 37:49) are correlated with each other across the ten sites, using a principal components analysis. Provide an appropriate plot and summary of the PCA results. What are the major patterns of variation across these sites?
# we can ask “what is the dominant axis of variation among these two variables?”. PCA answers this question by finding the the line through the data that explains the most variation in the data.
enviro.small = subset(zp, select = c(37:49)) #env variables
enviro.small = na.omit(enviro.small)
pca = princomp(enviro.small, cor = TRUE) #each variable has the same total variance
summary(pca)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 2.4591664 1.5418092 1.3316315 1.09389399 0.90423622
## Proportion of Variance 0.4651923 0.1828597 0.1364033 0.09204647 0.06289563
## Cumulative Proportion 0.4651923 0.6480519 0.7844552 0.87650164 0.93939726
## Comp.6 Comp.7 Comp.8 Comp.9
## Standard deviation 0.69506314 0.49771057 0.198369951 0.132877250
## Proportion of Variance 0.03716252 0.01905506 0.003026972 0.001358182
## Cumulative Proportion 0.97655978 0.99561485 0.998641818 1.000000000
## Comp.10 Comp.11 Comp.12 Comp.13
## Standard deviation 2.499503e-08 1.077036e-08 0 0
## Proportion of Variance 4.805782e-17 8.923119e-18 0 0
## Cumulative Proportion 1.000000e+00 1.000000e+00 1 1Info on the axes
Loadings are vectors that tell us the direction of the principal components
# What are these axes? We can see with loadings():
loadings(pca)##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mangrove.shoreline 0.329 0.352 0.270 0.122 0.205
## non.mangrove.shoreline 0.197 -0.517 0.163 0.192 -0.155 0.150
## percent.mangrove.shoreline 0.238 0.420 0.266 0.316 0.226
## pond.area 0.351 0.235 0.219 0.289 0.376
## d13C -0.112 -0.465 -0.193 0.498 0.509 -0.352
## C -0.213 0.315 0.236 0.478 0.133 -0.629
## d15N 0.141 0.300 0.281 -0.581 -0.344 0.181 -0.173
## N 0.574 -0.142 0.471 -0.578
## relief -0.354 0.121 0.151 -0.228 0.381 0.290
## percent.trees -0.340 0.231 -0.179 0.269 0.194 0.382
## precip -0.357 0.399 0.143
## percent.impervious 0.330 0.179 -0.314 0.312
## percent.developed 0.334 0.167 -0.304 0.313
## Comp.8 Comp.9 Comp.10 Comp.11 Comp.12 Comp.13
## mangrove.shoreline 0.406 0.658 0.176
## non.mangrove.shoreline -0.415 -0.237 -0.379 0.456
## percent.mangrove.shoreline 0.362 -0.176 -0.179 -0.348 -0.416 0.238
## pond.area -0.351 -0.368 0.138 0.369 -0.343
## d13C 0.206 0.103 0.179
## C -0.224 -0.195 0.213 -0.120
## d15N 0.309 -0.221 0.202 0.243 0.196
## N -0.213 -0.119 0.108
## relief -0.552 -0.135 0.217 -0.373 -0.200
## percent.trees -0.111 0.650 -0.111 0.304
## precip 0.536 -0.335 -0.367 0.336 0.150
## percent.impervious -0.338 0.429 -0.294 -0.244 -0.458
## percent.developed -0.347 -0.197 0.138 -0.161 0.509 0.438
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.077 0.077 0.077 0.077 0.077 0.077 0.077 0.077 0.077
## Cumulative Var 0.077 0.154 0.231 0.308 0.385 0.462 0.538 0.615 0.692
## Comp.10 Comp.11 Comp.12 Comp.13
## SS loadings 1.000 1.000 1.000 1.000
## Proportion Var 0.077 0.077 0.077 0.077
## Cumulative Var 0.769 0.846 0.923 1.000Biplots
The idea of the biplot is to summarize PCA results by plotting two things at once, using two PC axes: 1) the scores of the observations, 2) the relationship between the original variables and the PC axes.
# distance biplot, scaled so distance bt the observation scores approximates the euclidean distance bt the observations in raw data
biplot(pca, pch = 21, col = c('grey', 'blue'), scale = 1) 
# correlation biplot, scald so angle bt the vectors approximates the correlation bt those vectors, and angle bt vector and axis approximates the correlation bt the vector and PC component
biplot(pca, pch = 21, col = c('grey', 'blue'), scale = 0)
The positive principle component 1 shows percent mangrove shoreline, while percent trees and precipitation are on the negative end. For the positive principle component 2, percent trees, d15N, and carbon are associated with mangrove shoreline.
Unconstrained ordination
- Next, characterize the major axes of variation in zooplankton community composition (i.e., an unconstrained ordination). Consider how to transform and standardize the community data before the analysis. Provide an appropriate plot and interpret what it means. Finally, use envfit() to assess which environmental variables are correlated with the major axes of community composition.
PCOA
Patterns of abundance from the species data
community <- zp[, 8:35]
community = wisconsin(sqrt(community))
pcoa = dbrda(community ~ 1, dist = "bray", add = TRUE)
summary(pcoa)##
## Call:
## dbrda(formula = community ~ 1, distance = "bray", add = TRUE)
##
## Partitioning of Lingoes adjusted squared Bray distance:
## Inertia Proportion
## Total 10.15 1
## Unconstrained 10.15 1
##
## Eigenvalues, and their contribution to the Lingoes adjusted squared Bray distance
##
## Importance of components:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7
## Eigenvalue 1.2350 1.0707 0.83962 0.7934 0.70387 0.63825 0.57181
## Proportion Explained 0.1217 0.1055 0.08276 0.0782 0.06938 0.06291 0.05636
## Cumulative Proportion 0.1217 0.2273 0.31003 0.3882 0.45761 0.52052 0.57689
## MDS8 MDS9 MDS10 MDS11 MDS12 MDS13 MDS14
## Eigenvalue 0.54642 0.43864 0.38778 0.34680 0.30952 0.26152 0.24439
## Proportion Explained 0.05386 0.04324 0.03822 0.03418 0.03051 0.02578 0.02409
## Cumulative Proportion 0.63075 0.67398 0.71221 0.74639 0.77690 0.80268 0.82677
## MDS15 MDS16 MDS17 MDS18 MDS19 MDS20 MDS21
## Eigenvalue 0.22801 0.20307 0.19340 0.17104 0.16247 0.1481 0.14193
## Proportion Explained 0.02247 0.02002 0.01906 0.01686 0.01601 0.0146 0.01399
## Cumulative Proportion 0.84924 0.86926 0.88832 0.90518 0.92119 0.9358 0.94979
## MDS22 MDS23 MDS24 MDS25 MDS26
## Eigenvalue 0.1238 0.11084 0.10613 0.096128 0.072561
## Proportion Explained 0.0122 0.01093 0.01046 0.009475 0.007152
## Cumulative Proportion 0.9620 0.97291 0.98337 0.992848 1.000000NMDS
ord = metaMDS(community, dist = "bray", trymax = 20)## Run 0 stress 0.2638881
## Run 1 stress 0.2709061
## Run 2 stress 0.2670613
## Run 3 stress 0.2623633
## ... New best solution
## ... Procrustes: rmse 0.07993697 max resid 0.3252058
## Run 4 stress 0.2603333
## ... New best solution
## ... Procrustes: rmse 0.09585003 max resid 0.2922797
## Run 5 stress 0.279046
## Run 6 stress 0.2862657
## Run 7 stress 0.2694452
## Run 8 stress 0.2612046
## Run 9 stress 0.25988
## ... New best solution
## ... Procrustes: rmse 0.07839082 max resid 0.2405313
## Run 10 stress 0.3027268
## Run 11 stress 0.2698684
## Run 12 stress 0.2676447
## Run 13 stress 0.2918601
## Run 14 stress 0.2833772
## Run 15 stress 0.2715812
## Run 16 stress 0.2839804
## Run 17 stress 0.2761451
## Run 18 stress 0.2711986
## Run 19 stress 0.2841238
## Run 20 stress 0.2731261
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 19: stress ratio > sratmaxOrdination plot
plot(ord, type = 'text')
NMDS with groups
ord = metaMDS(community, distance = "bray", k = 2, trymax = 100)## Run 0 stress 0.2638881
## Run 1 stress 0.3225712
## Run 2 stress 0.2754784
## Run 3 stress 0.2768309
## Run 4 stress 0.2693917
## Run 5 stress 0.2934097
## Run 6 stress 0.2762871
## Run 7 stress 0.2735869
## Run 8 stress 0.3023357
## Run 9 stress 0.2639279
## ... Procrustes: rmse 0.1243625 max resid 0.3713615
## Run 10 stress 0.261243
## ... New best solution
## ... Procrustes: rmse 0.07313193 max resid 0.2087073
## Run 11 stress 0.2927377
## Run 12 stress 0.2655343
## Run 13 stress 0.2689849
## Run 14 stress 0.269774
## Run 15 stress 0.2636711
## Run 16 stress 0.2647576
## Run 17 stress 0.2811615
## Run 18 stress 0.2705186
## Run 19 stress 0.2681594
## Run 20 stress 0.2621046
## Run 21 stress 0.2639279
## Run 22 stress 0.2914551
## Run 23 stress 0.292336
## Run 24 stress 0.2648466
## Run 25 stress 0.2657895
## Run 26 stress 0.2796564
## Run 27 stress 0.2705833
## Run 28 stress 0.2994401
## Run 29 stress 0.2788736
## Run 30 stress 0.2603331
## ... New best solution
## ... Procrustes: rmse 0.07748421 max resid 0.2377563
## Run 31 stress 0.2681278
## Run 32 stress 0.2750737
## Run 33 stress 0.2683304
## Run 34 stress 0.2636711
## Run 35 stress 0.2713779
## Run 36 stress 0.2825172
## Run 37 stress 0.2900956
## Run 38 stress 0.2711214
## Run 39 stress 0.2697739
## Run 40 stress 0.2849141
## Run 41 stress 0.2781967
## Run 42 stress 0.2758213
## Run 43 stress 0.2690763
## Run 44 stress 0.2769613
## Run 45 stress 0.2688503
## Run 46 stress 0.2811175
## Run 47 stress 0.2896262
## Run 48 stress 0.2712281
## Run 49 stress 0.3167436
## Run 50 stress 0.274421
## Run 51 stress 0.2701198
## Run 52 stress 0.2681594
## Run 53 stress 0.2685896
## Run 54 stress 0.2696823
## Run 55 stress 0.2688699
## Run 56 stress 0.3008795
## Run 57 stress 0.2871188
## Run 58 stress 0.2648923
## Run 59 stress 0.2679956
## Run 60 stress 0.2878619
## Run 61 stress 0.2811244
## Run 62 stress 0.2670618
## Run 63 stress 0.2774689
## Run 64 stress 0.2916758
## Run 65 stress 0.2743585
## Run 66 stress 0.2865504
## Run 67 stress 0.2637298
## Run 68 stress 0.2709611
## Run 69 stress 0.2825425
## Run 70 stress 0.2817292
## Run 71 stress 0.3097083
## Run 72 stress 0.2924929
## Run 73 stress 0.2747979
## Run 74 stress 0.2754819
## Run 75 stress 0.3129392
## Run 76 stress 0.2794459
## Run 77 stress 0.2693317
## Run 78 stress 0.2748493
## Run 79 stress 0.2956015
## Run 80 stress 0.2936993
## Run 81 stress 0.2694449
## Run 82 stress 0.2969432
## Run 83 stress 0.3123909
## Run 84 stress 0.2966427
## Run 85 stress 0.2826387
## Run 86 stress 0.27289
## Run 87 stress 0.2784042
## Run 88 stress 0.276714
## Run 89 stress 0.2681693
## Run 90 stress 0.3108758
## Run 91 stress 0.2670612
## Run 92 stress 0.2760528
## Run 93 stress 0.2612934
## Run 94 stress 0.2825401
## Run 95 stress 0.2822517
## Run 96 stress 0.2696767
## Run 97 stress 0.2654942
## Run 98 stress 0.2759008
## Run 99 stress 0.2999009
## Run 100 stress 0.2658851
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 99: stress ratio > sratmaxord2 = metaMDS(community, distance = "bray", k = 3, trymax = 100)## Run 0 stress 0.1807737
## Run 1 stress 0.178888
## ... New best solution
## ... Procrustes: rmse 0.08776531 max resid 0.2813128
## Run 2 stress 0.1807976
## Run 3 stress 0.179379
## ... Procrustes: rmse 0.1226783 max resid 0.3960824
## Run 4 stress 0.1807977
## Run 5 stress 0.1744734
## ... New best solution
## ... Procrustes: rmse 0.04129465 max resid 0.1716513
## Run 6 stress 0.1827903
## Run 7 stress 0.1793785
## Run 8 stress 0.1840165
## Run 9 stress 0.1807971
## Run 10 stress 0.1829529
## Run 11 stress 0.1803415
## Run 12 stress 0.1803418
## Run 13 stress 0.1810127
## Run 14 stress 0.1744732
## ... New best solution
## ... Procrustes: rmse 0.0002513952 max resid 0.0007045174
## ... Similar to previous best
## Run 15 stress 0.1744733
## ... Procrustes: rmse 0.0001709031 max resid 0.0004122433
## ... Similar to previous best
## Run 16 stress 0.1744738
## ... Procrustes: rmse 0.0007799854 max resid 0.001764863
## ... Similar to previous best
## Run 17 stress 0.1793786
## Run 18 stress 0.1970492
## Run 19 stress 0.1744737
## ... Procrustes: rmse 0.0008018415 max resid 0.002194503
## ... Similar to previous best
## Run 20 stress 0.1803415
## *** Best solution repeated 4 timesord3 = metaMDS(community, distance = "bray", k = 4, trymax = 100)## Run 0 stress 0.1248321
## Run 1 stress 0.1248322
## ... Procrustes: rmse 0.0008161029 max resid 0.002761999
## ... Similar to previous best
## Run 2 stress 0.1279573
## Run 3 stress 0.1279587
## Run 4 stress 0.1248326
## ... Procrustes: rmse 0.0003532379 max resid 0.001205388
## ... Similar to previous best
## Run 5 stress 0.1328362
## Run 6 stress 0.1248426
## ... Procrustes: rmse 0.002203773 max resid 0.006483054
## ... Similar to previous best
## Run 7 stress 0.1351379
## Run 8 stress 0.1279575
## Run 9 stress 0.1248326
## ... Procrustes: rmse 0.0005490005 max resid 0.001723035
## ... Similar to previous best
## Run 10 stress 0.1279592
## Run 11 stress 0.1248323
## ... Procrustes: rmse 0.0001097922 max resid 0.0003357115
## ... Similar to previous best
## Run 12 stress 0.1248325
## ... Procrustes: rmse 0.0005073617 max resid 0.001886701
## ... Similar to previous best
## Run 13 stress 0.1331088
## Run 14 stress 0.1248319
## ... New best solution
## ... Procrustes: rmse 0.0002225003 max resid 0.000671714
## ... Similar to previous best
## Run 15 stress 0.1248322
## ... Procrustes: rmse 0.0005358971 max resid 0.001980788
## ... Similar to previous best
## Run 16 stress 0.1248321
## ... Procrustes: rmse 0.0004808528 max resid 0.001740836
## ... Similar to previous best
## Run 17 stress 0.1353847
## Run 18 stress 0.127958
## Run 19 stress 0.1279662
## Run 20 stress 0.1328632
## *** Best solution repeated 3 timesenvfit()
form = as.formula(paste("ord ~", paste(names(enviro.small),collapse = "+" )))
fit = envfit(form, data = zp, na.rm = T)
ordiplot(ord, display = "sites", type = 'n', main = 'nmds 2D')
points(ord, col = 'grey')
plot(fit, col = 'red', arrow.mul = 0.7) 
#fit = envfit(ord4 ~ nitrate + phosphate + silicate + PAR + chl + temp, data = enviro, na.rm = T) ordiplot(ord4, display = "sites", type = 'n', main = 'nmds 4D')
#fitted(ord.taxa)[, 1] fitted(ord.taxa)[, 2] points(ord4, col = 'grey')
#plot(fit, col = 'red', arrow.mul = 0.7)Constrained ordination
- Although envfit() is a useful tool, it does not directly ask how well community composition can be explained by a set of environmental variables. Use a constrained ordination to ask this question. How much of variation in community composition can be explained by the set of environmental variables? How much can be explained by the first two axes? What do those axes represent (visualizing with a triplot will help)?
## PCOA
community = wisconsin(sqrt(community))
pcoa = dbrda(community ~ 1, dist = "bray", add = TRUE)
# dbRDA
formdbRDA = as.formula(paste("community ~", paste(names(enviro.small),collapse = "+" )))
channel.dbrda = dbrda(formdbRDA, data = zp, dist = "bray", add = TRUE) ##
## Some constraints or conditions were aliased because they were redundant. This
## can happen if terms are constant or linearly dependent (collinear):
## 'percent.trees', 'precip', 'percent.impervious', 'percent.developed'plot(channel.dbrda)
PERMANOVA
- The authors were interested in whether zooplankton composition is affected by the presence of mangroves, and whether composition is affected by the presence of a fishpond. Use a PERMANOVA to test these questions. Considering all the results together, what have you learned about this system?
zpcomp <- zp %>%
group_by(Site, Mangrove, Fishpond) %>%
summarise(across(8:35, mean))## `summarise()` has grouped output by 'Site', 'Mangrove'. You can override using
## the `.groups` argument.Interaction
# permanova.interaction = adonis2(butter.wisc ~ Treatment*Year + Block, data = butterwide, by = 'margin')
# butter.counts = butterwide[,7:ncol(butterwide)]
# butter.wisc = wisconsin(sqrt(butter.counts))
community.counts = zpcomp[, 4:ncol(zpcomp)]
community.wisc = wisconsin(sqrt(community.counts))
permanova.interaction = adonis2(community.wisc ~ Mangrove * Fishpond, data = zpcomp, by = 'margin')
permanova.interaction## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = community.wisc ~ Mangrove * Fishpond, data = zpcomp, by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## Mangrove:Fishpond 1 0.13631 0.05124 1.1021 0.337
## Residual 16 1.97886 0.74382
## Total 19 2.66041 1.00000No interaction
permanova.no.interaction = adonis2(community.wisc ~ Mangrove + Fishpond, data = zpcomp, by = 'margin')
permanova.no.interaction## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = community.wisc ~ Mangrove + Fishpond, data = zpcomp, by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## Mangrove 1 0.14321 0.05383 1.1510 0.323
## Fishpond 1 0.40202 0.15111 3.2311 0.001 ***
## Residual 17 2.11517 0.79506
## Total 19 2.66041 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The mangrove:fishpond interaction and presence of a fishpond impacts zooplankton composition, showing us that this system is impacted by different environmental variables.