class: inverse, middle, center, the-title-slide, title-slide .title[ # Workshop in R ] .subtitle[ ## Washington State University 2024 ] .author[ ### Rob Smith ] .date[ ### 2024-04-04 ] --- class: inverse middle center big-subsection # Good afternoon! --- class: inverse middle center big-subsection # Overview --- # Overview 1\. Introductions 2\. Basic syntax and structures 3\. Usage and examples --- class: inverse middle center big-subsection # Introductions --- # Software installation 1\. Download and install R from: https://cran.r-project.org/ 2\. Optional: download and install Rstudio from: https://posit.co/download/rstudio-desktop/ 3\. Install packages (only once!) ```r pkg <- c('vegan', 'labdsv', 'ade4', 'ecodist', 'ape', 'picante') install.packages(pkg) ``` 4\. Load packages ```r sapply(pkg, require, character.only = TRUE) ``` 5\. Load some utility functions ```r source("https://github.com/phytomosaic/esa2022/raw/main/R/utils_color.R") ``` --- # Get the data ```r u <- 'https://github.com/phytomosaic/esa2023/raw/main/data/veg.rda' u <- url(u) load(u) ``` --- # Everything in R is one of: 1\. object 2\. function 3\. environment -- # Kinds of objects in R 1\. NULL 2\. scalar 3\. vector 4\. matrix or array 5\. data.frame 6\. list --- # Basic R syntax ```r ### assignment x <- 10 # objects: a scalar m <- matrix(c(1:8,NA), nrow=3, ncol=3) # objects: a 3x3 matrix f <- function(a, b) { a + b } # functions e <- new.env() # environments ### object structure str(x) # describe structure head(x) # peek at top few rows ### indexing values in an object m[1,1] # matrix row 1, column 1 lst$the_name # list item by name (data.frame column) ``` --- # Types and classes ```r ### examples x1 <- 1 x2 <- 1.0 x3 <- '1' x4 <- as.factor(x3) x5 <- TRUE x6 <- x5 * 99 ### examine classes list(x1,x2,x3,x4,x5,x6) lapply(list(x1,x2,x3,x4,x5,x6), class) ``` --- # Handling missing elements ```r ### create matrix object m <- matrix(c(1:8,NA), nrow=3, ncol=3) # a 3x3 matrix ### identify missings anyNA(m) is.na(m[1,1]) is.na(m[2,2]) is.na(m[3,3]) is.na(m) ### replacing missings m[is.na(m)] <- 777 # assign one value to any missing values ``` --- # Getting data into R Ideally, everything is "text" ```r ### read csv x <- read.csv('/path/to/your/file.csv') ### load binary RDA from file load('/path/to/your/file.rda') ### load binary data from package data('thedata') ``` --- # The Matrix  --- # The data ### Mafragh, Algeria vegetation data We will use one core dataset throughout this tutorial. The `veg` dataset gives information about spatial coordinates, species, environment, traits, and phylogeny for plants on the Mafragh coastal plain in North Africa. Specifically, `veg` is a list containing five items: - `xy` 97 observations of 2 spatial coordinates - `spe` 97 observations of 56 plant species - `env` 97 observations of 11 soil environmental variables - `tra` 12 traits for 56 plant species - `phy` phylogeny for 56 plant species --- # Load data ```r load('./data/veg.rda') # load the object (contains multiple objects) xy <- veg$xy # spatial spe <- veg$spe # species env <- veg$env # environment tra <- veg$tra # traits phy <- veg$phy # phylogeny rm(veg) # cleanup ls() # list objects now in this local environment ``` ``` ## [1] "colvec" "env" "get_palette" "phy" "pkg" ## [6] "snk" "spe" "tra" "xy" ``` -- <br></br> #### How would you examine what's inside one of these objects? --- class: inverse middle center big-subsection # Transformations --- # Transformations ```r ### Species: express abundances as 0/1 presence/absence # (spe > 0) * 1 ### Species: express abundances on log10 scale spe <- data.frame(log10(spe + 1)) ### Environment: express variables in [0,1] range env <- data.frame(vegan::decostand(env, 'range')) ### Traits: express traits in [0,1] range tra <- data.frame(vegan::decostand(tra, 'range')) ``` -- #### Why do we transform community data? --- class: inverse middle center big-subsection # Outliers --- # Outliers #### Define multivariate outlier function ```r `outliers` <- function (x, mult=2, method='bray') { d <- as.matrix(vegan::vegdist(x, method=method, binary=F, diag=T, upper=T)) diag(d) <- as.numeric(1) # avoid zero-multiplication m <- apply(d, 2, mean) # site means z <- scale(m) # z-scores data.frame(mean_dist = m, z = z, is_outlier = abs(z) >= mult) } ``` -- ```r o <- outliers(spe, mult=2) head(o, 3) ``` ``` ## mean_dist z is_outlier ## 1 0.8487643 0.3107496 FALSE ## 2 0.8283680 -0.2028542 FALSE ## 3 0.8214233 -0.3777314 FALSE ``` ```r which(o$is_outlier) ``` ``` ## [1] 32 45 95 97 ``` --- # Test validity of species matrix Species matrix must have no missings, no empty SUs, no empty species. ```r !anyNA(spe) # expect TRUE, no missing values all(rowSums(spe, na.rm=T) > 0) # expect TRUE, no empty sites all(colSums(spe, na.rm=T) > 0) # expect TRUE, no empty species ``` --- # Visualize data .pull-left.w50[ ```r ### spatial *plot(xy, pch=19, col='grey') ### species vegan::tabasco(spe, col=get_palette()) ### environment vegan::tabasco(env, col=get_palette()) ### traits vegan::tabasco(tra, col=get_palette()) ### phylogeny plot(phy, cex=0.6, no.margin=TRUE) ``` ] .pull-right.w50[ <!-- --> ] --- # Visualize data .pull-left.w50[ ```r ### spatial plot(xy, pch=19, col='grey') ### species *vegan::tabasco(spe, col=get_palette()) ### environment vegan::tabasco(env, col=get_palette()) ### traits vegan::tabasco(tra, col=get_palette()) ### phylogeny plot(phy, cex=0.6, no.margin=TRUE) ``` ] .pull-right.w50[ <!-- --> ] --- # Visualize data .pull-left.w50[ ```r ### spatial plot(xy, pch=19, col='grey') ### species vegan::tabasco(spe, col=get_palette()) ### environment *vegan::tabasco(env, col=get_palette()) ### traits vegan::tabasco(tra, col=get_palette()) ### phylogeny plot(phy, cex=0.6, no.margin=TRUE) ``` ] .pull-right.w50[ <!-- --> ] --- # Visualize data .pull-left.w50[ ```r ### spatial plot(xy, pch=19, col='grey') ### species vegan::tabasco(spe, col=get_palette()) ### environment vegan::tabasco(env, col=get_palette()) ### traits *vegan::tabasco(tra, col=get_palette()) ### phylogeny plot(phy, cex=0.6, no.margin=TRUE) ``` ] .pull-right.w50[ <!-- --> ] --- # Visualize data .pull-left.w50[ ```r ### spatial plot(xy, pch=19, col='grey') ### species vegan::tabasco(spe, col=get_palette()) ### environment vegan::tabasco(env, col=get_palette()) ### traits vegan::tabasco(tra, col=get_palette()) ### phylogeny *plot(phy, cex=0.6, no.margin=TRUE) ``` ] .pull-right.w50[ <!-- --> ] --- class: inverse middle center big-subsection # Diversity --- # Diversity measures #### How do you define "diversity"? ```r ### Gamma (regional) diversity (gamma <- sum(colSums(spe) > 0)) ``` ``` ## [1] 56 ``` ```r ### Alpha (per-site) diversity alpha <- rowSums(spe > 0) # within-site (avgalpha <- mean(alpha)) # average within-site ``` ``` ## [1] 6.268041 ``` ```r ### Beta (among-site) diversity: Whittaker's (beta <- gamma / avgalpha - 1) ``` ``` ## [1] 7.934211 ``` --- # Diversity measures (&#x3B2;-diversity) ```r ### 1 -- proportion of zeros in the matrix (independent of abundance) propzero <- sum(spe < .Machine$double.eps) / prod(dim(spe)) cat('Proportion of zeros in matrix:', propzero, '\n') ``` ``` ## Proportion of zeros in matrix: 0.8880707 ``` ```r ### 2 -- "dust bunny index" of McCune and Root (2015) (uses abundances) dbi <- 1 - mean(as.matrix(vegan::decostand(spe, method='max'))) cat('Dust bunny index:', dbi, '\n') ``` ``` ## Dust bunny index: 0.9202329 ``` ```r ### 3 -- pairs of SUs that don't share species z <- vegan::no.shared(spe) propnoshare <- sum(z) / length(z) cat('', propnoshare, 'proportion of site-pairs share no species in common\n') ``` ``` ## 0.441366 proportion of site-pairs share no species in common ``` --- class: inverse middle center big-subsection # Dissimilarities --- # Dissimilarities .pull-left.w48[ #### Species ```r d <- vegdist(spe, method='bray', binary=T) tabasco(as.matrix(d), col=get_palette()) ``` <!-- --> ] .pull-right.w48[ #### Environment ```r E <- vegdist(env, method='euc', binary=F) tabasco(as.matrix(E), col=get_palette()) ``` <!-- --> ] --- # Loss of sensitivity problem .pull-left.w60[ How dissimilar are SUs A and C? They share no species in common... ``` ## sp1 sp2 sp3 sp4 sp5 ## suA 1 1 0 0 0 ## suB 0 1 1 1 0 ## suC 0 0 0 1 1 ``` ] -- .pull-left.w60[ ...Bray-Curtis will max out at 1.0... ``` ## suA suB ## suB 0.6 ## suC 1.0 0.6 ``` ] -- .pull-left.w60[ ...but `stepacross()` replaces "too long" distances with shortest path. ``` ## suA suB ## suB 0.6 ## suC 1.2 0.6 ``` ] --- # Stepacross adjustment <img src="index_files/figure-html/double-zero-05-1.svg" style="display: block; margin: auto;" /> --- # Stepacross adjustment ```r D <- stepacross(d, 'shortest', toolong = 1) plot(d, D, xlab = 'Original', ylab = 'Stepacross') ``` <img src="index_files/figure-html/diss-02-invisible-eval-1.svg" style="display: block; margin: auto;" /> --- class: inverse middle center big-subsection # Ordination --- # Ordination: unconstrained ### Three different algorithms ```r m1 <- metaMDS(D, k=2, maxit=250, try=100, trymax=101, trace=0) # NMS m2 <- cmdscale(D, k=2, add=T) # PCoA m3 <- prcomp(spe) # PCA ``` <!-- --> --- # dbRDA: the Swiss army knife ```r ### dissimilarities D_euc <- vegdist(spe, 'euc') # Euclidean D_chi <- vegdist(spe, 'chi') # Chi-sq D_bc <- vegdist(spe, 'bray') # Bray-Curtis (Sorensen) ### PCA: PCO based on Euclidean distances (see `stats::prcomp`) pca <- dbrda(D_euc ~ 1) ### CA/RA: PCO based on Chi-sq distances (see `vegan::cca`) ca <- dbrda(D_chi ~ 1) ### PCO: generalizes to any dissimilarity (see `labdsv::pco`) pco <- dbrda(D_bc ~ 1) ### RDA: constrained form of PCA (see `vegan::rda`) rda <- dbrda(D_euc ~ k2o + mg, data = env) ### CCA: constrained form of CA (see `vegan::cca`) cca <- dbrda(D_chi ~ k2o + mg, data = env) ### dbRDA: constrained form of PCO; can handle neg eigenvalues dbr <- dbrda(D_bc ~ k2o + mg, data = env, add = 'lingoes') ``` --- # dbRDA: the Swiss army knife .pull-left.w85[ <img src="index_files/figure-html/ord-swissarmy-02-1.svg" style="display: block; margin: auto;" /> ] .pull-right.w15[ ##### What are the points? ##### What is the ordination space? ##### What does distance between points mean? ] --- # Scaling scores Using `scores(x, scaling = ...)`: - `scaling = 1` &mdash; focus on sites, scale site scores by `\(\lambda_i\)` - `scaling = 2` &mdash; focus on species, scale species scores by `\(\lambda_i\)` - `scaling = 3` &mdash; **symmetric scaling**, scale both scores by `\(\sqrt{\lambda_i}\)` - `scaling = -1` &mdash; as above, but for `rda()` get correlation scores - `scaling = -2` &mdash; for `cca()` multiply results by `\(\sqrt{(1/(1-\lambda_i))}\)` - `scaling = -3` &mdash; this is Hill's scaling - `scaling < 0` &mdash; for `rda()` divide species scores by species' `\(\sigma\)` - `scaling = 0` &mdash; raw scores ...where `\(\lambda_i\)` is the *i*th eigenvalue. <br></br> ###### Credit: Gavin Simpson --- # NMS: a sensible default ```r nms <- vegan::metaMDS(D, k = 2, maxit = 500, try = 500, trymax = 501) plot(nms) # barebones -- what would you do to interpret this? stressplot(nms) ``` <img src="index_files/figure-html/ord-nms-hide-1.svg" style="display: block; margin: auto;" /> --- class: inverse middle center big-subsection # Group clustering --- # Group clustering ### Hierarchical: Ward's clustering ```r k <- 7 # specify number of groups cl <- hclust(D, method='ward.D2') # clustering solution grp <- cutree(cl, k) # group memberships plot(cl) ; rect.hclust(cl, k, border=rainbow(k)) # plot the dendrogram ``` <img src="index_files/figure-html/cluster-01-1.svg" style="display: block; margin: auto;" /> --- # Group clustering ### Non-hierarchical: fuzzy clustering ```r install.packages(vegclust) require(vegclust) k <- 7 # specify number of groups cl <- vegclust::vegclustdist(D, mobileMemb=k) # clustering solution grp_f <- cl$memb # fuzzy membership grp_c <- vegclust::defuzzify(cl, 'cut', alpha=0.8) # crisp membership ``` --- class: inverse middle center big-subsection # Group differences --- # Group differences #### First, define and visualize groups (arbitrary example) ```r plot(m1$points, pch=NA, asp=1) # visualize on NMS text(m1$points, labels=grp, col=as.numeric(grp)) # group memberships ordispider(m1, groups=grp, col=1:k) # group centroids ``` <img src="index_files/figure-html/permanova-grp-plot-1.svg" style="display: block; margin: auto;" /> --- # Group differences #### PERMANOVA: test for differences in multivariate *centroid* ```r adonis2(D ~ grp, permu=99) ``` ``` ## Permutation test for adonis under reduced model ## Terms added sequentially (first to last) ## Permutation: free ## Number of permutations: 99 ## ## adonis2(formula = D ~ grp, permutations = 99) ## Df SumOfSqs R2 F Pr(>F) ## grp 1 9.510 0.20354 24.278 0.01 ** ## Residual 95 37.214 0.79646 ## Total 96 46.724 1.00000 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` -- #### Recall how we defined the groups in `grp`; are you comfortable with this hypothesis test? --- # Group differences #### PERMDISP: test for differences in multivariate *dispersion* ```r permutest(betadisper(D, grp), pairwise=TRUE, permu=99) ``` ``` ## Df Sum Sq Mean Sq F N.Perm Pr(>F) ## Groups 6 0.4708883 0.07848139 4.364745 99 0.01 ## Residuals 90 1.6182676 0.01798075 NA NA NA ``` ``` ## Permuted p-values: ``` ``` ## 1-2 1-3 1-4 1-5 1-6 1-7 2-3 2-4 2-5 2-6 2-7 3-4 3-5 3-6 3-7 4-5 ## 0.62 0.06 0.04 0.02 0.04 0.30 0.45 0.19 0.07 0.46 0.37 0.01 0.01 0.47 0.02 0.43 ## 4-6 4-7 5-6 5-7 6-7 ## 0.05 0.43 0.02 0.16 0.10 ``` --- # PERMANOVA extensions #### Can use categorical and/or continuous predictors ```r adonis2(D ~ grp + elevation + clay + silt + conductivity, permu=99) ``` #### Blocked design: permutations must occur *within* strata ```r blk <- factor(letters[sample(rep(1:12,len=nrow(env)))]) # arbitrary 'blocks' pm <- how(nperm=999) # setup permutation object setBlocks(pm) <- blk # permute only *within* blocks adonis2(D ~ grp, permu=pm) # correct test ``` --- class: inverse middle center big-subsection # Group indicator species --- # Group indicator species #### Real groups ```r iv <- labdsv::indval(spe, grp) # indicator species analysis for *real* groups summary(iv) # IndVal observed ``` ``` ## cluster indicator_value probability ## halimione_portulacoides 1 0.1854 0.047 ## atriplex_prostrata 1 0.2731 0.009 ## carlina_racemosa 1 0.4144 0.002 ## centaurium_spicatum 1 0.4172 0.002 ## scilla_autumnalis 1 0.5288 0.001 ## hordeum_marinum 1 0.6991 0.001 ## lolium_rigidum 1 0.7014 0.001 ## alisma_plantago_aquatica 2 0.3194 0.020 ## schoenoplectus_litoralis 2 0.7283 0.001 ## aeluropus_littoralis 3 0.3929 0.008 ## bolboschoenus_maritimus 3 0.5295 0.001 ## phalaris_coerulescens 4 0.2832 0.037 ## beta_vulgaris 4 0.4280 0.005 ## carlina_lanata 5 0.2305 0.035 ## plantago_coronopus 5 0.2854 0.009 ## trifolium_squamosum 5 0.2966 0.023 ## medicago_intertexta 5 0.3037 0.005 ## arthrocnemum_macrostachyum 5 0.3666 0.013 ## juncus_maritimus 6 0.3400 0.005 ## sarcocornia_fruticosa 6 0.7541 0.001 ## damasonium_alisma 6 0.7925 0.001 ## rumex_bucephalophorus 7 0.2344 0.045 ## stachys_marrubiifolia 7 0.2819 0.013 ## convolvulus_tricolor 7 0.3333 0.012 ## picris_echioides 7 0.4019 0.012 ## convolvulus_arvensis 7 0.4085 0.002 ## sonchus_oleraceus 7 0.4731 0.004 ## avena_fatua 7 0.6111 0.001 ## borago_officinalis 7 0.6111 0.001 ## diplotaxis_erucoides 7 0.7222 0.001 ``` --- # Group indicator species #### Random groups ```r rnd <- sample(grp, length(grp), replace=T) # define random groups by bootstrapping ivr <- labdsv::indval(spe, rnd) # indicator species analysis for *random* groups summary(ivr) # IndVal expected at random ``` ``` ## cluster indicator_value probability ## drimia_numidica 6 0.2903 0.008 ``` -- #### Null expectation, setting alpha = 0.05 ```r ceiling(ncol(spe) * 0.05) ``` ``` ## [1] 3 ``` --- class: inverse middle center big-subsection # Community traits --- # Community traits ```r head(tra, 4) ``` ``` ## anemogamous autogamous entomogamous annual biennial ## arisarum_vulgare 0 0 1 0 0 ## alisma_plantago_aquatica 0 0 1 0 0 ## damasonium_alisma 0 0 1 1 1 ## asphodelus_aetivus 0 0 1 0 0 ## perennial lfp min_height max_height bfp ## arisarum_vulgare 1 0.6444426 0.4114081 0.1879018 0.0000000 ## alisma_plantago_aquatica 1 0.4065980 0.4114081 0.4362945 0.6666667 ## damasonium_alisma 1 0.6444426 0.2342243 0.1099155 0.5555556 ## asphodelus_aetivus 1 0.5374933 1.0000000 0.5462100 0.2777778 ## spikiness hairy_leaves ## arisarum_vulgare 0 0 ## alisma_plantago_aquatica 0 0 ## damasonium_alisma 0 0 ## asphodelus_aetivus 0 0 ``` --- # Community traits .pull-left.w50[ ```r ### Euclidean trait dissimilarity ### (traits already scaled 0-1) Dt <- dist(tra, method='euc') ### calculate trait SES of mean ### pairwise distances in sites ses <- picante::ses.mpd(spe, Dt, null.model='richness') ### plot SES across a nutrient gradient u <- ifelse(ses$mpd.obs.p < 0.05,2,1) plot(ses$mpd.obs.z ~ env$k, ylab='Trait SES(MPD)', xlab='Soil potassium', col=u) abline(lm(ses$mpd.obs.z ~ env$k)) # regression abline(h=0, lty=2) # random-traits line text(0.9, 1, 'Divergent') text(0.9, -4, 'Convergent') ``` ] .pull-right.w50[ <!-- --> ] --- # Community weighted means .pull-left.w50[ #### (Weighted) mean trait value per SU ```r ### function to make CWM matrix `makecwm` <- function (spe, tra) { spe <- as.matrix(spe) tra <- as.matrix(tra) `stdz` <- function(x) { (x - min(x, na.rm=TRUE)) / diff(range(x, na.rm=TRUE)) } tra <- apply(tra, MARGIN = 2, FUN = stdz) awt <- spe %*% tra # abund-weighted totals awt / rowSums(spe, na.rm=TRUE) # CWM matrix } ### make the CWM traits matrix cwm <- data.frame(makecwm(spe, tra)) ### visualize tabasco(cwm, col=get_palette()) ``` ] .pull-right.w50[ <!-- --> ] --- # Ordination of CWM traits .pull-left.w50[ #### Admits nonlinear trait-enviro relationship ```r ### NMS based on abundance-weighted traits m <- metaMDS(cwm, 'altGower', k=2, trace=0) ### SUs sized relative to a leaf trait plot(m, cex=cwm$lfp) ### overlay enviro variable o <- ordisurf(m, env$k, col=2, add=T) ``` ##### What are the points? ##### What is the ordination space? ##### What does distance between points indicate? ] .pull-right.w50[ <!-- --> ] --- # Fourth-corner analysis and RLQ ### the RLQ method ```r require(ade4) o_spe <- dudi.coa(spe, F) # Correspondence Analysis o_env <- dudi.hillsmith(env, F, row.w = o_spe$lw) o_tra <- dudi.hillsmith(tra, F, row.w = o_spe$cw) r <- rlq(o_env, o_spe, o_tra, F) plot(r) summary(r) randtest(r) ``` --- # Fourth-corner analysis and RLQ <img src="index_files/figure-html/rlq-02-1.svg" style="display: block; margin: auto;" /> --- # Fourth-corner analysis and RLQ #### Trait correlations ```r fourthcorner.rlq(r, type='Q.axes') ``` ``` ## Fourth-corner Statistics ## ------------------------ ## Permutation method Comb. 2 and 4 ( 999 permutations) ## ## Adjustment method for multiple comparisons: holm ## call: fourthcorner.rlq(xtest = r, typetest = "Q.axes") ## ## --- ## ## Test Stat Obs Std.Obs Alter Pvalue ## 1 AxcR1 / anemogamous r 0.3305888147 5.69292462 two-sided 0.001 ## 2 AxcR2 / anemogamous r -0.0458582587 -0.77960338 two-sided 0.444 ## 3 AxcR1 / autogamous r -0.2245440975 -2.08793024 two-sided 0.03 ## 4 AxcR2 / autogamous r 0.0516421890 1.11307304 two-sided 0.259 ## 5 AxcR1 / entomogamous r -0.3305888147 -5.69292462 two-sided 0.001 ## 6 AxcR2 / entomogamous r 0.0458582587 0.77960338 two-sided 0.444 ## 7 AxcR1 / annual r -0.2710550381 -2.52925916 two-sided 0.007 ## 8 AxcR2 / annual r -0.0984989149 -1.67108085 two-sided 0.089 ## 9 AxcR1 / biennial r -0.1276065206 -1.13982521 two-sided 0.261 ## 10 AxcR2 / biennial r -0.0641373702 -1.67709458 two-sided 0.092 ## 11 AxcR1 / perennial r 0.2705213859 2.48889593 two-sided 0.007 ## 12 AxcR2 / perennial r 0.0340545959 0.52989683 two-sided 0.59 ## 13 AxcR1 / lfp r 0.0737622617 1.75201478 two-sided 0.086 ## 14 AxcR2 / lfp r 0.0756252046 1.26281835 two-sided 0.206 ## 15 AxcR1 / min_height r -0.0002066631 -0.03086926 two-sided 0.976 ## 16 AxcR2 / min_height r 0.1004326492 1.63046869 two-sided 0.109 ## 17 AxcR1 / max_height r 0.0051630195 0.07343021 two-sided 0.951 ## 18 AxcR2 / max_height r 0.0544587825 1.29079312 two-sided 0.201 ## 19 AxcR1 / bfp r 0.2324526902 2.12347283 two-sided 0.026 ## 20 AxcR2 / bfp r -0.0909484362 -1.50781687 two-sided 0.145 ## 21 AxcR1 / spikiness r 0.1779342409 1.75271427 two-sided 0.081 ## 22 AxcR2 / spikiness r 0.0636401007 1.40307023 two-sided 0.167 ## 23 AxcR1 / hairy_leaves r -0.1842685566 -1.69495302 two-sided 0.091 ## 24 AxcR2 / hairy_leaves r -0.0097666513 -0.20355949 two-sided 0.842 ## Pvalue.adj ## 1 0.024 * ## 2 1 ## 3 0.57 ## 4 1 ## 5 0.024 * ## 6 1 ## 7 0.154 ## 8 1 ## 9 1 ## 10 1 ## 11 0.154 ## 12 1 ## 13 1 ## 14 1 ## 15 1 ## 16 1 ## 17 1 ## 18 1 ## 19 0.52 ## 20 1 ## 21 1 ## 22 1 ## 23 1 ## 24 1 ## ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- # Fourth-corner analysis and RLQ #### Environment correlations ```r fourthcorner.rlq(r, type='R.axes') ``` ``` ## Fourth-corner Statistics ## ------------------------ ## Permutation method Comb. 2 and 4 ( 999 permutations) ## ## Adjustment method for multiple comparisons: holm ## call: fourthcorner.rlq(xtest = r, typetest = "R.axes") ## ## --- ## ## Test Stat Obs Std.Obs Alter Pvalue ## 1 clay / AxcQ1 r 0.158920982 2.44189264 two-sided 0.014 ## 2 silt / AxcQ1 r -0.061153914 -0.96433329 two-sided 0.326 ## 3 sand / AxcQ1 r -0.115820713 -1.79748542 two-sided 0.077 ## 4 k2o / AxcQ1 r 0.268313148 4.11998443 two-sided 0.001 ## 5 mg / AxcQ1 r 0.160253315 2.39412055 two-sided 0.014 ## 6 na_100g / AxcQ1 r 0.353355511 5.34613443 two-sided 0.001 ## 7 k / AxcQ1 r 0.449746890 6.69092568 two-sided 0.001 ## 8 conductivity / AxcQ1 r 0.316892896 4.83476275 two-sided 0.001 ## 9 retention / AxcQ1 r 0.118177718 1.72785058 two-sided 0.094 ## 10 na_l / AxcQ1 r 0.280641685 4.24166475 two-sided 0.001 ## 11 elevation / AxcQ1 r -0.262101010 -2.67217392 two-sided 0.002 ## 12 clay / AxcQ2 r 0.079573076 1.48111204 two-sided 0.145 ## 13 silt / AxcQ2 r -0.155231156 -2.92952837 two-sided 0.006 ## 14 sand / AxcQ2 r 0.037032254 0.71340595 two-sided 0.472 ## 15 k2o / AxcQ2 r 0.026201484 0.50296940 two-sided 0.602 ## 16 mg / AxcQ2 r -0.076592437 -1.55057344 two-sided 0.113 ## 17 na_100g / AxcQ2 r -0.032223166 -0.65956317 two-sided 0.522 ## 18 k / AxcQ2 r 0.108332685 0.91226907 two-sided 0.386 ## 19 conductivity / AxcQ2 r -0.059251811 -1.19077435 two-sided 0.245 ## 20 retention / AxcQ2 r -0.060635579 -1.20040885 two-sided 0.225 ## 21 na_l / AxcQ2 r -0.082406453 -1.67655037 two-sided 0.093 ## 22 elevation / AxcQ2 r 0.003330203 0.03305252 two-sided 0.968 ## Pvalue.adj ## 1 0.21 ## 2 1 ## 3 0.924 ## 4 0.022 * ## 5 0.21 ## 6 0.022 * ## 7 0.022 * ## 8 0.022 * ## 9 1 ## 10 0.022 * ## 11 0.034 * ## 12 1 ## 13 0.096 . ## 14 1 ## 15 1 ## 16 1 ## 17 1 ## 18 1 ## 19 1 ## 20 1 ## 21 1 ## 22 1 ## ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- class: inverse middle center big-subsection # Community phylogenetics --- # Community phylogenetics: basic plotting ```r par(mfrow=c(1,3), mar=c(0,0,0,0), oma=c(0,0,0,0)) # plotting parameters plot(phy, cex=0.75, no.margin=T) # basic plotting plot(phy, cex=0.75, tip.color=colvec(tra$lfp), no.margin=T) # color labels by traits plot(phy, cex=0.75, label.offset=48, no.margin=T) # symbolize trait values at tips tiplabels(pch=21, bg = c(tra$annual), adj=0) tiplabels(pch=21, bg = c(tra$biennial), adj=15) tiplabels(pch=21, bg = c(tra$perennial), adj=30) ``` <img src="index_files/figure-html/phyloplot-01-1.svg" style="display: block; margin: auto;" /> --- # Phylogenetic diversity - alpha .pull-left.w50[ ```r ### phylogenetic distance matrix Dp <- cophenetic(phy) ### Faith's phylogenetic diversity: ### total branch length connecting ### all species per site fpd <- pd(spe, phy) ### mean pairwise distance mpd <- ses.mpd( spe, Dp, null.model='independentswap') ### mean nearest taxon distance mnd <- ses.mntd( spe, Dp, null.model='independentswap') ### bring all together phy_div <- cbind(fpd, mpd = mpd$mpd.obs.z, mntd = mnd$mntd.obs.z) ``` ] .pull-right.w47[ ``` ## PD SR mpd mntd ## 1 767.32 6 0.92152467 1.08223033 ## 2 551.60 4 1.00127789 0.89956948 ## 3 364.80 3 -0.74186332 -0.52641802 ## 4 694.38 6 -0.23269719 0.08822563 ## 5 395.80 3 -0.26034867 0.33585163 ## 6 323.60 3 -1.03913225 -1.66716963 ## 7 485.76 4 -0.17221377 -0.18773460 ## 8 932.32 9 0.04080037 0.08182749 ## 9 476.80 4 -1.01997248 -0.09066436 ## 10 579.60 6 -1.65459527 -1.13524065 ## 11 879.94 10 -0.48729928 -1.97117056 ## 12 817.38 7 -0.09923170 0.53220249 ``` ] --- # Phylogenetic diversity - beta ```r ### correlation between phylogenetic and taxonomic beta-diversity Dp <- picante::phylosor(spe, phy) # phylogenetic distances protest(Dp, D) # procrustes correlation ``` ``` ## ## Call: ## protest(X = Dp, Y = D) ## ## Procrustes Sum of Squares (m12 squared): 0.2736 ## Correlation in a symmetric Procrustes rotation: 0.8523 ## Significance: 0.001 ## ## Permutation: free ## Number of permutations: 999 ``` --- # Phylogenetic signal ```r sapply(tra, FUN=function(j){ names(j) <- rownames(tra) round(picante::phylosignal(j, ape::multi2di(phy)), 4)}) ``` ``` ## anemogamous autogamous entomogamous annual biennial ## K 3.2136 0.3401 3.2136 0.3466 0.4564 ## PIC.variance.obs 6e-04 0.0039 6e-04 0.0057 0.0016 ## PIC.variance.rnd.mean 0.0066 0.0047 0.0066 0.0069 0.0026 ## PIC.variance.P 0.001 0.251 0.001 0.148 0.158 ## PIC.variance.Z -5.3551 -0.6838 -5.316 -1.1289 -1.0962 ## perennial lfp min_height max_height bfp spikiness ## K 0.3951 0.3512 0.3451 0.3043 0.4064 0.9436 ## PIC.variance.obs 0.005 0.0014 9e-04 9e-04 9e-04 0.001 ## PIC.variance.rnd.mean 0.0069 0.0017 0.0011 0.001 0.0012 0.0033 ## PIC.variance.P 0.04 0.17 0.266 0.464 0.124 0.001 ## PIC.variance.Z -1.7276 -0.8997 -0.6953 -0.2408 -1.083 -2.6172 ## hairy_leaves ## K 0.5273 ## PIC.variance.obs 0.0024 ## PIC.variance.rnd.mean 0.0044 ## PIC.variance.P 0.007 ## PIC.variance.Z -2.3532 ``` --- class: inverse middle center big-subsection # Community spatial analysis --- # Mantel test ```r E <- dist(xy) # spatial distance matrix vegan::mantel(D, E, method = 'spearman') # spearman *rank* correlation ``` ``` ## ## Mantel statistic based on Spearman's rank correlation rho ## ## Call: ## vegan::mantel(xdis = D, ydis = E, method = "spearman") ## ## Mantel statistic r: 0.4093 ## Significance: 0.001 ## ## Upper quantiles of permutations (null model): ## 90% 95% 97.5% 99% ## 0.0257 0.0350 0.0426 0.0495 ## Permutation: free ## Number of permutations: 999 ``` ```r ecodist::mantel(D ~ E, mrank = T) # alternative with useful bootstrap CIs ``` --- # Mantel correlogram ```r plot(vegan::mantel.correlog(D, E, cutoff=F, r.type='spearman', nperm=99, mult='holm')) ``` <img src="index_files/figure-html/mantel-04-1.svg" style="display: block; margin: auto;" /> --- # Multiple regression on distance matrices .pull-left[ ```r # species dissimilarities are related # to space (extremely weakly) and # potassium (moderately) ecodist::MRM( D ~ dist(xy) + dist(env$k)) ``` ``` ## $coef ## D pval ## Int 0.627953130 1.000 ## dist(xy) 0.001755839 0.001 ## dist(env$k) 0.301671934 0.001 ## ## $r.squared ## R2 pval ## 0.1994291 0.0010000 ## ## $F.test ## F F.pval ## 579.5514 0.0010 ``` ] .pull-right[ ```r # abundance of cosmopolitan bulrush # is NOT related to space, but is # related to potassium (moderately) ecodist::MRM( dist(spe$bolboschoenus) ~ dist(xy) + dist(env$k)) ``` ``` ## $coef ## dist(spe$bolboschoenus) pval ## Int 0.201552720 1.000 ## dist(xy) 0.000223733 0.091 ## dist(env$k) 0.411223852 0.001 ## ## $r.squared ## R2 pval ## 0.06277826 0.00100000 ## ## $F.test ## F F.pval ## 155.8368 0.0010 ``` ] --- class: inverse middle center big-subsection --- class: inverse middle center big-subsection # Wrap up --- # Acknowledgments ### ESA Vegetation Section Martin Dovciak Carissa Brown Morgan D. Frost Anita Thompson ### Further details https://ecol.shinyapps.io/esa_tutorial/ --- class: inverse middle center big-subsection # More?? --- # PCNM PCNM: principle coordinates of neighbor matrices ```r E <- dist(xy) # euclidean distances between sites pc <- pcnm(E) # principal coordinates of neighbor matrices ``` <!-- --> --- # PCNM ```r rs <- rowSums(spe) / sum(spe) # sites weighted by abundances pc <- pcnm(E, w=rs) # *weighted* PCNMs vp <- varpart(D, env, vegan::scores(pc)) # variation partitioning ``` <img src="index_files/figure-html/pcnm-03b-1.svg" style="display: block; margin: auto;" /> --- # PCNM .pull-left[ ```r ### broad and local structure pc_broad <- vegan::scores(pc)[,1:10] pc_local <- vegan::scores(pc)[,11:59] vp <- varpart(D, env, pc_broad, pc_local) ``` ] .pull-right[ <img src="index_files/figure-html/pcnm-04b-1.svg" style="display: block; margin: auto;" /> ] --- # PCNM .pull-left[ ```r ### soil physical, soil chemistry, and spatial predictors vp <- varpart( D, # dissimilarities ~ clay + sand + silt, # soil fractions ~ mg + k + conductivity + na_l, # soil chemistry vegan::scores(pc), # spatial predictors data = env) # environmental dataset ``` ] .pull-right[ <img src="index_files/figure-html/pcnm-05b-1.svg" style="display: block; margin: auto;" /> ] --- class: inverse middle center big-subsection # Conclusion --- class: inverse middle center big-subsection ---