Differences

This shows you the differences between two versions of the page.

--- r_workshop10 [2015/09/08 14:56]
zofia.taranu [Workshop 10: Programming in R]
+++ r_workshop10 [2021/10/13 20:05] (current)
lsherin
@@ Line 1: / Line 1: @@
-======= QCBS R Workshops =======
+<WRAP group>
+<WRAP centeralign>
+<WRAP important>
+<wrap em> __MAJOR UPDATE__ </wrap>
-[[http://qcbs.ca/|{{:logo_text.png?nolink&500|}}]]
+<wrap em> As of Fall 2021, this wiki has been discontinued and is no longer being actively developed. </wrap>
-This series of [[r|10 workshops]] walks participants through the steps required to use R for a wide array of statistical analyses relevant to research in biology and ecology.
+<wrap em> All updated materials and announcements for the QCBS R Workshop Series are now housed on the [[https://r.qcbs.ca/workshops/r-workshop-10/|QCBS R Workshop website]]. Please update your bookmarks accordingly to avoid outdated material and/or broken links. </wrap>
-These open-access workshops were created by members of the QCBS both for members of the QCBS and the larger community.
-====== Workshop 10: Programming in R ======
-Developed by: Johanna Bradie, Sylvain Christin, Ben Haller, Guillaume Larocque
+<wrap em> Thank you for your understanding, </wrap>
-**Summary:** This workshop focuses on basic programming in R.  In this workshop, you will learn how to use control flow (for loops, if, while) methods to prevent code repetition, facilitate organization and run simulations.  In addition, you will learn to write your own functions, and tips to program efficiently. The last part of the workshop will discuss packages that were not covered in this workshop series, but that may be of interest to participants.
+<wrap em> Your QCBS R Workshop Coordinators. </wrap>
-Link to associated Prezi: [[https://prezi.com/xuu2rphp5wg4/qcbs-r-workshop-8/|Prezi]]
+</WRAP>
+</WRAP>
+<WRAP clear></WRAP>
-Download the R script for this lesson:
+======= QCBS R Workshops =======
-  * [[http://qcbs.ca/wiki/_media/qcbs_workshop10_r_code.r| R Script]]  *
+[[http://qcbs.ca/|{{:logo_text.png?nolink&500|}}]]
+This series of [[r|10 workshops]] walks participants through the steps required to use R for a wide array of statistical analyses relevant to research in biology and ecology.
+These open-access workshops were created by members of the QCBS both for members of the QCBS and the larger community.
-===== Learning Objectives =====
+//The content of this workshop has been peer-reviewed by several QCBS members. If you would like to suggest modifications, please contact the current series coordinators, listed on the main wiki page//
-  - Flow control
+<wrap em>IMPORTANT NOTICE: MAJOR UPDATES</wrap>
-  - Writing functions in R
-  - Speeding up your code
-  - Useful packages for biologists
-===== 1. Flow Control =====
+**March 2021 update:** This wiki has been discontinued and is no longer being actively developed or updated. Updated materials for the QCBS R Workshop Series are now housed on the QCBS R Workshop [[https://github.com/QCBSRworkshops/workshop10|GitHub page]].
-Flow control allows you to run the same series of commands multiple times and subject to specified conditions.  In this section, you will learn how to:
+Available materials include;
+  - The [[https://qcbsrworkshops.github.io/workshop10/pres-en/workshop10-pres-en.html|Rmarkdown presentation]] for this workshop;
+  - An [[https://qcbsrworkshops.github.io/workshop10/book-en/workshop10-script-en.R|R script]] that follows the presentation;
+  - [[https://qcbsrworkshops.github.io/workshop10/book-en/index.html|Written materials]] that accompany the presentation in bookdown format.
-  - Execute statements **conditionally** using: if, if/else
-  - Execute statements **multiple times** using: for loops, while loops, repeat loops
-  - **Modify loop** execution using: break statements, next statements
-==== if and if/else statements ====
+====== Workshop 10: Advanced multivariate analyses ======
-if and if/else statements are good for:
+Developed by: Monica Granados, Emmanuelle Chrétien, Bérenger Bourgeois, Amanda Winegardner and Xavier Giroux-Bougard. (Material in R script adapted from: Borcard, Gillet & Legendre (2011). //Numerical Ecology with R//. Springer New York.)
-  * checking for problems or violations of assumptions
-  * treating different rows of your data frame differently
-  * testing for the existence of a file or variable
-**Syntax**
+**Summary:** In this workshop, you will learn how to perform advanced multivariate analyses on community data. This workshop concentrates on constrained methods such as redundancy analysis (RDA), multivariate regression tree (MRT) and linear discriminant analysis (LDA) to explore how environmental variables may be driving patterns in species assemblage across sites.
-<code rsplus>
+**Link to new [[https://qcbsrworkshops.github.io/workshop10/workshop10-en/workshop10-en.html|Rmarkdown presentation]]**
-if (condition) {
-  expression      # The expression can be any command that you would like R to perform.
-}
-if (condition) {
+Link to old [[https://prezi.com/zzqqe4gcq80g/|Prezi presentation]]
-  expression
-} else {
-  expression
-}
-</code>
+Download the R script, packages and data required for this workshop:
+  *  [[http://qcbs.ca/wiki/_media/multivar2_e.r| R Script]]
+  *  [[http://qcbs.ca/wiki/_media/DoubsEnv.csv|DoubsEnv data]]
+  *  [[http://qcbs.ca/wiki/_media/DoubsSpe.csv|DoubsSpe data]]
+  *  [[http://qcbs.ca/wiki/_media/classifyme.csv|Test data for linear discriminant analysis]]
+  *  [[http://qcbs.ca/wiki/_media/mvpart_1.6-2.tar.gz|mvpart package]]
+  *  [[http://qcbs.ca/wiki/_media/MVPARTwrap_0.1-9.tar.gz|MVPARTwrap package]]
+  *  [[http://qcbs.ca/wiki/_media/rdaTest_1.10.tar.gz|rdaTest package]]
-Throughout this workshop, we will be making use of **logical operators**.
-Note that to test whether something is equal to something else, you must use a "==".
-| == | equal to |
+Make sure to load the following packages (see how in the R script):
-| != | not equal to |
+  *  [[http://cran.r-project.org/web/packages/vegan/index.html|vegan (for multivariate analyses)]]
-| !x | not x |
+  *  [[http://cran.r-project.org/web/packages/vegan/index.html|labdsv (for identification of significant indicator species in the multivariate regression tree analysis)]]
-| < | less than |
+  *  [[http://cran.r-project.org/web/packages/vegan/index.html|plyr (classification for linear discriminant analysis)]]
-| < = | less than or equal to |
+  *  [[http://cran.r-project.org/web/packages/vegan/index.html|MASS (for linear discriminant analysis)]]
-| > | greater than |
+  *  mvpart
-| >= | greater than or equal to |
+  *  MVPARTwrap
-| x & y | x AND y |
+  *  rdatest
-| x%%|%%y| x OR y |
-| isTRUE(x) | test if X is true |
-For example,
-<code rsplus>
+<code rsplus | Load the required packages and functions>
-if ((2+2) == 4) {
+install.packages("vegan")
-  print("Arithmetic works.")
+install.packages("mvpart")
-}
+install.packages("labdsv")
+install.packages("plyr")
+install.packages("MASS")
-if ((2+1) == 4) {
+# For the two following packages, upload the file provided on the wiki page.
-  print("Arithmetic works.")
+# To do so, go to Packages tab on the bottom right panel of R Studio
-}
+# Click on Install Packages
-</code>
+# Choose to install from Package Archive file and upload these two files
+install.packages("MVPARTwrap")
+install.packages("rdaTest")
-Curly brackets { } are used so that R knows to expect more input.  When using brackets, R waits to evaluate the command until the brackets have been closed.  If the curly brackets are not used, R may not behave as you are expecting.  For example, try:
+library(vegan)
+library(mvpart)
-<code rsplus>
+library(MVPARTwrap)
-if ((2 + 1) == 4) print("Arithmetic works.")
+library(rdaTest)
-else print("Houston, we have a problem.")
+library(labdsv)
+library(plyr)
+library(MASS)
 </code>
-The else statement doesn't work because R evaluates the first line without knowing your command is incomplete.
-Instead, use:
-<code rsplus>
-if ((2 + 2) == 4) {
-  print("Arithmetic works.")  # R does not evaluate this expression yet because the bracket isn't closed.
-} else {
-  print("Houston, we have a problem.")
-}  # Since all brackets are now closed, R will evaluate the commands.
-</code>
-Note that if and if/else test a single condition.  If you want to test a vector of conditions (and get a vector of results), you can use the ifelse function:
+======Introduction to advanced multivariate analyses======
-For example,
+The previous workshop presented the basics of multivariate analyses:
-<code rsplus>
+    * how to choose appropriate distance metrics and transformations
-a <- 1:10
+    * hierarchical clustering
-ifelse(a > 5, "yes", "no")
+    * unconstrained ordinations
-</code>
+    * Principal component analysis
+    * Principal coordinate Analysis
+    * Correspondence analysis
+    * Nonmetric multidimensional scaling
-You can also use the ifelse function within another function to carry out an operation only under certain conditions:
+The present workshop builds on this knowledge, and will focus on constrained analyses.  All the methods overviewed during the introductory workshop allowed to find patterns in the community composition data or in the descriptors, but not to explore how environmental variables could be driving these patterns. With constrained analyses, such as redundancy analysis (RDA), linear discriminant analysis (LDA) and multivariate regression tree (MRT), one can describe and predict relationships between community composition data and environmental variables.
-For example,
-<code rsplus>
-a <- (-4):5
-sqrt(ifelse(a >= 0, a, NA))
-</code>
+======Getting started with data======
------
+We will continue to use the Doubs river datasets for this workshop. “DoubsSpe.csv” is a data frame of fish community data where the first column contains site names from 1 to 30 and the remaining columns are fish taxa. The taxa columns are populated by fish abundance data (counts). “DoubsEnv.csv” is a data frame of environmental data for the same sites contained in the fish community data frame. Again, the first column contains site names from 1 to 30. The remaining columns contain measurements for 11 abiotic variables. Note that data used in ordination analyses is generally in [[http://en.wikipedia.org/wiki/Wide_and_narrow_data|wide-format]].
-**Exercise 1**
-<code rsplus>
+<code rsplus | Load the Doubs Species and Environmental Data>
-Paws <- "cat"
+#Species community data frame (fish abundance): “DoubsSpe.csv”
-Scruffy <- "dog"
+spe<- read.csv(file.choose(), row.names=1)
-Sassy <- "cat"
+spe<- spe[-8,] #Site number 8 contains no species and so row 8 (site 8) is removed. Be careful to
-animals <- c(Paws, Scruffy, Sassy)
+#only run this command line once as you are overwriting "spe" each time.
-</code>
-. Use an if statement to print "meow" if Paws is a "cat".
-++++ Exercise 1.1 : Answer|
+#Environmental data frame: “DoubsEnv.csv”
-<code rsplus>
+env<- read.csv(file.choose(), row.names=1)
-if(Paws == "cat") {
+env<- env[-8,] #Remove corresponding abiotic data for site 8 (since removed from fish data).
-  print("meow")
+#Again, be careful to only run the last line once.
-}
 </code>
-++++
-. Use an if/else statement to print "woof" if you supply an object that is a "dog" and "meow" if it is not.  Try it out with Paws and Scruffy.
-++++ Exercise 1.2 : Answer|
+======1. Explore and prepare the data======
-<code rsplus>
-if(Scruffy == "dog") {
-  print("woof")
-} else {
-  print("meow")
-}
-if(Paws == "dog") {
+=====1.1 Species data =====
-  print("woof")
-} else {
-  print("meow")
-}
-</code>
-++++
-. Use the ifelse function to display "woof" for animals that are dogs and "meow" for animals that are cats.
+We can use summary functions to explore the “spe” data (fish community data) and discover things like the dimensions of the matrix, column headings and summary statistics for the columns. This is a review from Workshop 2.
-++++ Exercise 1.3 : Answer|
+<code rsplus | Explore DoubsSpe>
-<code rsplus>
+names(spe) #see names of columns in spe
-ifelse(animals == "dog", "woof", "meow")
+dim(spe) #dimensions of spe; number of columns and rows
+str(spe) #displays internal structure of objects
+head(spe) #first few rows of the data frame
+summary(spe) #summary statistics for each column; min value, median value, max value, mean value etc.
 </code>
-++++
------
-==== Loops ====
-Loops are good for:
+Look at the species’ distribution frequencies.
-  * doing something for every element of an object
-  * doing something until the processed data runs out
-  * doing something for every file in a folder
-  * doing something that can fail, until it succeeds
-  * iterating a calculation until it converges
-===for loops===
+<code rsplus | Species distribution of DoubsSpe data>
+#Species distribution
+(ab <- table(unlist(spe))) #note that when you put an entire line of code in brackets like this, the output for that operation is displayed right away in the R console
-The //for loop// is the most common type of loop.  Use a //for loop// to execute a block of code a known number of times.
+barplot(ab, las=1, xlab="Abundance class", ylab="Frequency", col=grey(5:0/5))
-**Syntax**
-<code rsplus>
-for (variable in sequence) {
-  expression
-}
 </code>
-Each time the series of commands in a loop are executed, it is known as an iteration.
+{{:spe_barplot.png?300|}}
+Can see that there is a high frequency of zeros in the abundance data.
-For example:
+See how many absences there are in the fish community data.
-<code rsplus>
+<code rsplus | Absences in fish data>
-for (i in 1:5) {
+sum(spe==0)
-  print(i)
-}
 </code>
-In this example, our sequence has 5 elements (1, 2, 3, 4, 5). R will then evaluate the expression 5 times. The variable i is specific to our loop and will be only accessible inside it. Its value will be each of the element of our sequence. Therefore, in the first iteration, R will replace each instance of i with 1.  In the second iteration i would be replaced with 2, and so on.
+Look at the proportion of zeros in the fish community data.
+<code rsplus | Proportion of zeros>
-The letter 'i' can be replaced with any variable name and the sequence can be almost anything, even a list of vectors.
+sum(spe==0)/(nrow(spe)*ncol(spe))
-Try:
-<code rsplus>
-for (m in 4:10) {
-  print(m * 2)
-}
-for (a in c("Hello", "R", "Programmers")) {
-  print(a)
-}
-for (z in 1:30) {
-  a <- rnorm(n = 1, mean = 5, sd = 2) # draw a value from a normal distribution with mean 5 and standard deviation 2
-  print(a)
-}
-elements <- list(1:3, 4:10)
-for (element in elements) {
-  print(element)
-}
 </code>
+The proportion of zeros in the dataset is ~0.5.
-Loops are often used to loop over a dataset.  We will use loops to perform functions on the CO2 dataset which is built in to R.  The CO2 dataset contains concentration and uptake values for plants located in Quebec and Mississippi that were exposed to a treatment ("chilled") or control ("nonchilled").  Note that this is the dataset we used in workshop 2.  The code below provides a couple of examples of how loops can be used.
+This is high but not uncommon for species abundance data. However, in order to avoid the use of double-zeros as indications of resemblance among sites, we will apply a transformation to the species data. Legendre and Gallagher (2001) proposed five pre-transformations of the species data, four of them being available in vegan in the function decostand().
-<code rsplus>
+The Hellinger transformation will be applied to the fish data. It expresses abundances as the square-root of their relative abundance at each site (Borcard et al. 2011).
-data(CO2)
-for (i in 1:length(CO2[,1])) { # for each row in the CO2 dataset
-  print(CO2$conc[i]) #print the CO2 concentration
-}
-for (i in 1:length(CO2[,1])) { # for each row in the CO2 dataset
+<code rsplus | Hellinger transformation>
-  if(CO2$Type[i] == "Quebec") { # if the type is "Quebec"
+spe.hel <- decostand(spe, method="hellinger") # you can also use method="hell"
-    print(CO2$conc[i]) #print the CO2 concentration }
-  }
-}
-# Tip 1 : to get the number of rows of a data frame, we can also use the function nrow
-for (i in 1:nrow(CO2)) { # for each row in the CO2 dataset
-  print(CO2$conc[i]) #print the CO2 concentration
-}
-# Tip 2 : If we want to perform operations on only the elements of one column, we can directly
-# iterate over it.
-for (i in CO2$conc) { # for every element of the concentration column of the CO2 dataset
-  print(i) # print the ith element
-}
 </code>
-The expression part of the loop can be almost anything and is usually a compound statement containing many commands.
-<code rsplus>
+=====1.2 Environmental data=====
-for (i in 4:5) { # for i in 4 to 5
+Explore the environmental data and create a panel of plots to compare collinearity:
-  print(colnames(CO2)[i])
-  print(mean(CO2[,i])) # print the mean of that column from the CO2 dataset
-}
-</code>
-Note that this could be done more quickly using apply(), but that wouldn't teach you about loops. We will talk about it later.
-===while loops and repeat loops ===
-//while loops// and //repeat loops// operate similarly to //for loops//.  Once you understand how //for loops// work, you should be able to use any type of loop.  You will see some examples of //while loops// and //repeat loops// in the next section.  Note that in many cases you can accomplish the same task many different ways-- either using for loops, while loops, or repeat loops as will be demonstrated below.
-=== Nested loops ===
-In some cases, you may want to use nested loops to accomplish a task.  When using nested loops, it is important to use different variables as counters for each of your loops (here we used i and n).
-<code rsplus>
-for (i in 1:5) {
-  for (n in 1:5) {
-    print (i*n)
-  }
-}
+<code rsplus | Explore Doubs Env data>
+names(env)
+dim(env)
+str(env)
+head(env)
+summary(env)
+pairs(env, main="Bivariate Plots of the Environmental Data" )
 </code>
------
+In this case, the environmental data (explanatory variables) are all in different units and need to be standardized prior to computing distance measures to perform most ordination analyses. Standardize the environmental data (11 variables) using the function decostand() in vegan.
-**Exercise 2**
-. You have realized that your tool for measuring uptake was not calibrated properly at Quebec sites and all measurements are 2 units higher than they should be.  Use a loop to correct these measurements for all Quebec sites.
-++++ Exercise 2 : Answer|
-<code rsplus>
-for (i in 1:length(CO2[,1])) {
-  if(CO2$Type[i] == "Quebec") {
-    CO2$uptake[i] <- CO2$uptake[i] - 2
-  }
-}
+<code rsplus | Decostand>
+env.z <- decostand(env, method="standardize")
+apply(env.z, 2, mean) # the data are now centered (means~0)
+apply(env.z, 2, sd)   # the data are now scaled (standard deviations=1)
 </code>
-++++
------
-Make sure you reload the data so that we are working with the raw data for the rest of the exercise:
-<file rsplus>
-data(CO2)
-</file>
-\\
-==== Loop Modifications ====
-Normally, loops iterate over and over until they finish.  To change this behavior, you can use **break** which breaks out of the loops execution entirely, or **next**, which stops executing the current iteration and jumps to the next iteration.
-For example,
+======2. Canonical analyses======
-<file rsplus>
+Described first by Rao (1964), canonical analysis is a generic term that for several types of statistical analyses sharing a common goal; to identify the relationship between a multivariate response table (matrix Y, generally describing the species composition of communities) and a multivariate explanatory table (matrix X, generally containing environmental descriptors) by combining ordination and regression concepts. Canonical analyses allow users to test ecological hypothesis concerning the environmental drivers of species composition. Among the diversity of canonical analysis, we will mainly focus  here on Redundancy Analysis (RDA).
-# Print the CO2 concentrations for "chilled" treatments and keep count of how many replications there were.
-count <- 0 # count is being set at zero so that we can edit this object in the loop to keep track of how many iterations were performed
+=====2.1 Redundancy analysis (RDA)=====
-for (i in 1:length(CO2[,1])) {
+Redundancy Analysis is a direct extension of multiple regression, as it models the effect of an explanatory matrix X (n x p) on a response matrix Y (n x m). This is done by preforming an ordination of Y to obtain ordination axes that are linear combinations of the variables in X. In RDA, ordination axes are calculating from a PCA of a matrix Yfit, computed by fitting the Y variables to X by multivariate linear regression. Note that the explanatory variables in X can be quantitative, qualitative or binary variables. Prior to RDA, explanatory variables in Y must be centered, standardized (if explanatory variables are not dimensionally homogeneous, i.e. in different units), transformed (to limit the skew of explanatory variables) or normalized (to linearize relationships) following the same principles as in PCA. Collinearity between the X variables should also be reduced before RDA.
-  if (CO2$Treatment[i] == "nonchilled") next #Skip to next iteration if treatment is nonchilled
-  count <- count + 1
-  print(CO2$conc[i])
-}
-print(count) # The count and print command were performed 42 times.
-# This could be equivalently written using a repeat loop:
+In order to obtain the best model of RDA, explanatory variables can be selected by forward, backward or stepwise selection that remove non-significant explanatory variables.
-count <- 0
+RDA involves two computational steps. In the first step, a matrix of fitted values Yfit is calculated trough the linear equation:
-i <- 0
+Yfit = X[X'X]-1 [X'Y]
-repeat {
-      i <- i + 1
-      if (CO2$Treatment[i] == "nonchilled") next  # skip this loop
-      count <- count + 1
-      print(CO2$conc[i])
-      if (i == length(CO2[,1])) break     # stop looping
-    }
-print(count)
-### This could also be written using a while loop:
+In the second step, a PCA of the fitted matrix Yfit is calculated (see equations in the figure below) and produce the canonical eigenvalues and eigenvectors together with the matrix Z containing the canonical axes. These canonical axes correspond to linear combinations of the explanatory variables in X. This linearity of the combinations of the X variables is a fundamental property of RDA. In the analysis of community composition, these canonical axes are interpreted as complex environmental gradients.
-i <- 0
+{{ :constrained_ord_diagram.png |}}
-count <- 0
-while (i < length(CO2[,1]))
-{
-  i <- i + 1
-  if (CO2$Treatment[i] == "nonchilled") next  # skip this loop
-  count <- count + 1
-  print(CO2$conc[i])
-}
-print(count)
-</file>
------
-**Exercise 3**
-. You have realized that your tool for measuring concentration didn't work properly.  At Mississippi sites, concentrations less than 300 were measured correctly but concentrations >= 300 were overestimated by 20 units.  Use a loop to correct these measurements for all Mississippi sites.
+Redundancy analysis as a two-step process (from Legendre and Legendre 2012)
+Several statistics can be computed from RDA:
-++++ Exercise 3 : Answer|
+- The R² measures the strength of the canonical relationship between Y and X by calculating the proportion of the variation of Y explained by the variables in X,
-<file rsplus>
+- The adjusted R² also measures the strength of the relationship between Y and X, but applies a correction of the R² to take into account the number of explanatory variables,
-for (i in 1:length(CO2[,1])) {
+- The F-statistic corresponds to an overall test of significance of an RDA by comparing the computed model to a null model. This test is based on the null hypothesis that the strength of the linear relationship calculated by the R² is not larger than the value that would be obtained for unrelated Y and X matrices of the same size. Note that F-statistics can also be used to sequentially test the significance of each canonical axis.
-  if(CO2$Type[i] == "Mississippi") {
-    if(CO2$conc[i] < 300) next
-    CO2$conc[i] <- CO2$conc[i] - 20
-  }
-}
-# Note : We could also have written it that way, which is more concise and clear
+In R, RDA can be computed using the function rda from package vegan, as follows:
-for (i in 1:nrow(CO2)) {
-  if(CO2$Type[i] == "Mississippi" && CO2$conc[i] >= 300) {
-    CO2$conc[i] <- CO2$conc[i] - 20
-  }
-}
-</file>
-++++
------
-Make sure you reload the data so that we are working with the raw data for the rest of the exercise:
+<code rsplus | rda() in vegan for RDA>
-<file rsplus>
+#Preparing the data prior to RDA
-data(CO2)
+env.z <- subset(env.z, select = -das) # remove the "distance from the source" variable
-</file>
+#Running the RDA
+?rda
+spe.rda <- rda(spe.hel~., data=env.z)
+### Extract the results
+summary(spe.rda, display=NULL)
+#The results are called using summary:
+summary(spe.rda, display=NULL) #display = NULL optional
+</code>
-**Using flow control to make a complex plot**
-The idea here is that we have a dataset we want to plot, with concentration and uptake values, but each point has a type (Quebec or Mississippi) and a treatment ("chilled" or "nonchilled") and we want to plot the points differently for these cases.
+The summary of the output looks like this:
-You can read more about mathematical typesetting with ?plotmath, and more about the way that different colors, sizes, rotations, etc. are used in ?par.
-<file rsplus>
+{{ :rda_output_1.png |}}
-head(CO2) # Look at the dataset
+These results contain the proportion of variance of Y explained by the X variables (constrained proportion, 73.41% here), the unexplained variance of Y (unconstrained proportion, 26.59% here) and then summarize the eigenvalues, the proportions explained and the cumulative proportion of each canonical axis (each canonical axis = each constraining variable, in this case, the environmental variables from env).
-unique(CO2$Type)
-unique(CO2$Treatment)
-# plot the dataset, showing each type and treatment as a different colour
+To select the significant explanatory variables, you can then perform a forward selection (or backwards or stepwise), using the ordiR2step() function (or using the forward.sel function of package packfor):
-plot(x=CO2$conc, y=CO2$uptake, type="n", cex.lab=1.4, xlab="CO2 concentration", ylab="CO2 uptake") # Type "n" tells R to not actually plot the points.
-for (i in 1:length(CO2[,1])) {
-  if (CO2$Type[i] == "Quebec" & CO2$Treatment[i] == "nonchilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="red",type="p")
-  }
-  if (CO2$Type[i] == "Quebec" & CO2$Treatment[i] == "chilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="blue")
-  }
-  if (CO2$Type[i] == "Mississippi" & CO2$Treatment[i] == "nonchilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="orange")
-  }
-  if (CO2$Type[i] == "Mississippi" & CO2$Treatment[i] == "chilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="green")
-  }
-}
-</file>
+<code rsplus | ordiR2step() for forward selection>
+### Select the significant explanatory variables by forward selection
+?ordiR2step
+ordiR2step(rda(spe.hel~1, data=env.z), scope= formula(spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
+env.signif <- subset(env.z, select = c("alt", "oxy", "dbo"))
+</code>
-Note, that there are other ways to create a complex plot.  ggplot is a useful package for creating complex plots that was covered in workshop 4.
------
+The resulting output reads: rda(formula = spe.hel ~ alt + oxy + dbo, data = env.z) with the proportion of variation explained by the three constraining variables being 0.59. So in this case, only three variables are retained by the forward selection, i.e. alt, oxy and dbo. These three variables can be placed in a new data frame env.signif to perform the new RDA including only significant X variables:
-**Exercise 4**
-. Generate a plot of showing concentration versus uptake where each plant is shown using a different colour point.  Bonus points for doing it with nested loops!
-++++ Exercise 4 : Answer|
+<code rsplus | RDA with the two retained variables>
-<file rsplus>
+spe.rda.signif <- rda(spe.hel~., data=env.signif)
-plot(x=CO2$conc, y=CO2$uptake, type="n", cex.lab=1.4,xlab="CO2 concentration", ylab="CO2 uptake") # Type "n" tells R to not actually plot the points.
+summary(spe.rda.signif, display=NULL)
+</code>
-plants <- unique(CO2$Plant)
-for (i in 1:length(CO2[,1])){
+The explanatory variables (altitude, oxygen and biological oxygen demand) now explain 59% of the variance in Y (species).
-  for (p in 1:length(plants)) {
+The adjusted R² of this RDA is calculated using the function RsquareAdj:
-    if (CO2$Plant[i] == plants[p]) {
-      points(CO2$conc[i], CO2$uptake[i], col=p, type="p")
-    }
-  }
-}
-</file>
-++++
------
-===== 2. Writing Functions =====
-==== Why write functions? ====
+<code rsplus | adjusted R2>
+(R2adj <- RsquareAdj(spe.rda.signif)$adj.r.squared)
+#Here the strength of the relationship between X and Y corrected for the number of X variables is 0.54.
+</code>
-Much of the heavy lifting in R is done by functions. They are useful for:
-  * performing a task repeatedly, but configurably
-  * making your code more readable
-  * make your code easier to modify and maintain
-  * sharing code between different analyses
-  * sharing code with other people
-  * modifying R’s built-in functionality
+The significance of the model and of each canonical axis can be tested using the function anova (note this is different from retaining significant variables as was done with forward selection, now we're testing the significance of the RDA axes):
-But what exactly is a function? A function is essentially a black box that transforms data. It takes variables as entry - called arguments - use R code to process them and then can optionally give back a return value.
-{{::fonction_schema_eng2.png|200}}
+<code rsplus | anova.cca for testing the significance of axes>
+?anova.cca
-==== How to write functions? ====
+anova.cca(spe.rda.signif, step=1000)
+anova.cca(spe.rda.signif, step=1000, by="axis")
-Here is the basic syntax of a function:
+#In this case, the RDA model is highly significant (p=0.001) as well as all three canonical axes.
-<code rsplus>
-function_name <- function(argument1, argument2, ...) {
-  expression...  # What we want the function to do
-  return(value)  # Optional. If you want to access to the result of your function
-}
 </code>
-=== Arguments ===
-Arguments are the entry values of your function. They are the information your function needs to be able to perform correctly. A function can have between 0 and an infinity of arguments.
+To visualize the results of the RDA, triplots can be drawn using the plot(). Note that as in PCA, users can create scaling 1 and scaling 2 triplots. In scaling 1, distance among objects approximate their Euclidean distances while in scaling 2, angles between variables X and Y reflect their correlation. Thus, scaling 1 triplots can be used to interpret distances among objects and scaling 2 triplots to interpret the relationships between X and Y. To plot the scaling 1 triplots of the RDA, the following code can be used:
-From a technical standpoint, arguments are variables like any other so you can use them as such. The only difference is that they are available only inside your function. Their value will be determined at the moment your function will be called.
-Let's start with a really basic function we will call **print_number** that will take as argument a number and will print it.
-<code rsplus>
+<code rsplus | plot RDAs>
-print_number <- function(number) {
+#Quick plots scaling 1 and 2
-  print(number)
+windows()
-}
+plot(spe.rda.signif, scaling=1, main="Triplot RDA (scaling 1)")
-</code>
+windows()
+plot(spe.rda.signif, scaling=2, main="Triplot RDA (scaling 2)")
-Now to use it we call it like any other function.
+#Advanced plots scaling 1
-<code rsplus>
+windows()
-print_number(2)
+plot(spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
-print_number(231)
+points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(spe.rda.signif, display="species", choices=c(1), scaling=1),
+       scores(spe.rda.signif, display="species", choices=c(2), scaling=1),
+       col="black",length=0)
+text(scores(spe.rda.signif, display="species", choices=c(1), scaling=1),
+     scores(spe.rda.signif, display="species", choices=c(2), scaling=1),
+     labels=rownames(scores(spe.rda.signif, display="species", scaling=1)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(spe.rda.signif, display="bp", choices=c(1), scaling=1),
+       scores(spe.rda.signif, display="bp", choices=c(2), scaling=1),
+       col="red")
+text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
+     scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
+     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)),
+     col="red", cex=1)
+#Advanced plots scaling 2
+windows()
+plot(spe.rda.signif, scaling=2, main="Triplot RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=2),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2,
+       scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2,
+       col="black",length=0)
+text(scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2.1,
+     scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2.1,
+     labels=rownames(scores(spe.rda.signif, display="species", scaling=2)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(spe.rda.signif, display="bp", choices=c(1), scaling=2),
+       scores(spe.rda.signif, display="bp", choices=c(2), scaling=2),
+       col="red")
+text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=2)+0.05,
+     scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)+0.05,
+     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)),
+     col="red", cex=1)
 </code>
-We can use more than one arguments. For example let's create a function that take a number1, add it to a number2, multiply the result by a number3 and then prints the result.
-<code rsplus>
+And the final triplots would look like this:
-operations <- function(number1, number2, number3) {
-  result <- (number1 + number2) * number3
-  print(result)
-}
-operations(1, 2, 3)
+{{ :doubs_rda1.png |}}
-operations(17, 23, 2)
+{{ :doubs_rda2.png |}}
-</code>
-The expression part of our function - its body - can be virtually anything. It can be single R statements, loops, if/else conditions etc.
+**Challenge 1**: Run an RDA of the mite environmental variables constraining the mite species abundances. Use:
------
-**Exercise 5**
-Using what you learned previously on flow control, create a function **print_animal** that takes an animal as argument and gives the following results :
+<code rsplus | Load the mite data>
+#Load the mite species and environmental data from vegan package
+data(mite)
+mite.spe<-mite
+mite.spe.hel <- decostand(mite.spe, method="hellinger")
-<code rsplus>
+data(mite.env)
-Scruffy <- "dog"
-Paws <- "cat"
-print_animal(Scruffy)
-[1] "woof"
-print_animal(Paws)
-[1] "meow"
 </code>
-++++ Exercise 5 : Answer|
-<file rsplus| Exercise5 Answer>
-print_animal <- function(animal) {
-  if (animal == "dog") {
-    print("woof")
-  } else if (animal == "cat") {
-    print("meow")
-  }
-}
-</file>
-++++
------
-Arguments can also be optional and be provided with a default value. This is useful when using a function often with the same settings as it prevents the need from writing all the arguments all the time but still provides the flexibility to be able to change it if needed.
+What are the significant explanatory variables? How much of the variation in species data is explain by significant explanatory variables (i.e. those that were selected)? What are the significant axes? What group(s) of sites can you identify? What species are related to each group(s) of sites?
-<code rsplus>
-operations <- function(number1, number2, number3=3) {
-  result <- (number1 + number2) * number3
-  print(result)
-}
-operations(1, 2, 3) # becomes equivalent to
+**Challenge 1**: Solution
-operations(1, 2)
-operations(1, 2, 2) # we can still change the value of number3 if needed
+<hidden>
+Your code probably looks something like this:
+<code rsplus | RDA of mite data>
+#Initial RDA with ALL of the environmental data
+mite.spe.rda<-rda(mite.spe.hel~., data=mite.env)
+#Select significant environmental variables
+ordiR2step(rda(mite.spe.hel~1, data=mite.env),
+           scope= formula(mite.spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
+#Create a new dataframe with only the significant variables that you identified above
+mite.env.signif <- subset(mite.env,
+                          select = c("WatrCont", "Shrub", "Substrate", "Topo", "SubsDens"))
+#Re-run the RDA with the significant variables and look at the summary
+mite.spe.rda.signif=rda(mite.spe~., data=mite.env.signif)
+summary(mite.spe.rda.signif, display=NULL)
+#Find the R2 adjusted of the model with the retained environmental variables
+(R2adj <- RsquareAdj(mite.spe.rda.signif)$adj.r.squared)
+#Determine the significant canonical (constrained) axes)
+anova.cca(mite.spe.rda.signif, step=1000)
+anova.cca(mite.spe.rda.signif, step=1000, by="axis")
+#Plot the RDA
+windows()
+plot(mite.spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(mite.spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+text(scores(mite.spe.rda.signif, display="species", choices=c(1), scaling=1),
+     scores(mite.spe.rda.signif, display="species", choices=c(2), scaling=1),
+     labels=rownames(scores(mite.spe.rda.signif, display="species", scaling=1)),
+     col="grey", cex=0.8)
+arrows(0,0,
+      scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1),
+      scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1),
+      col="red")
+text(scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
+     scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
+     labels=rownames(scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)),
+     col="red", cex=1)
 </code>
-R also provides a special argument "**...**". It allows you to tell R that your function will accept an indefinite number of arguments. This is useful for two main things:
+Five explanatory variables are significant: WatrCont, Shrub, Substrate, Topo, and SubsDens. The proportion of the variance explained by these variables is 32.08% and the adjusted R² of this model is 19.20%. While the first canonical is significant (p<0.001), the overall model is however not significant (p=0.107). Three groups of sites appear on the triplot. One group of sites with high water content (top right) have high abundance of species 9 and 25. Species 6, 12, 19 and 26 are related to another group of sites with hummock topography (bottom left). The last group of sites (top left) have not characteristic species and heterogeneous environmental conditions. See the triplot below:
-  * Pass on arguments to another function used inside your function. This allows you to be able to use all the arguments of other functions without having to define them when creating yours. Let's create a function from our previous example where we plot the CO2 uptake based on the concentration. Here we will plot with 2 different colors depending on the region. Parameters to plot() and points() will be passed on via "...".
+{{ :mite_rda1.png |}}
-<code rsplus>
+</hidden>
-plot.CO2 <- function(CO2, ...) {
-  plot(x=CO2$conc, y=CO2$uptake, type="n", ...)  # We use ... to pass on arguments to plot().
-  for (i in 1:length(CO2[,1])){
-     if (CO2$Type[i] == "Quebec") {
-       points(CO2$conc[i], CO2$uptake[i], col="red", type="p", ...) # same for points
-     } else if (CO2$Type[i] == "Mississippi") {
-       points(CO2$conc[i], CO2$uptake[i], col="blue", type="p", ...) # same for points()
-     }
-  }
-}
-plot.CO2(CO2, cex.lab=1.4, xlab="CO2 concentration", ylab="CO2 uptake")
+=====2.2 Partial RDA=====
-plot.CO2(CO2, cex.lab=1.4, xlab="CO2 concentration", ylab="CO2 uptake", pch=20)
-</code>
+Partial RDA is a special case of RDA in which the response variables Y are related to explanatory variables X in the presence of additional explanatory variables, W, called covariables. As in partial linear regression, the linear effect of X variables on the Y variables are adjusted for the effects of the covariables W. For this, a RDA of the covariables W on the response variables Y is first performed. The residuals of this RDA are then extracted, i.e. a matrix Yres|W containing the Y response variables in which the effect of W were removed. The partial RDA finally correspond to the RDA of X on Yres|W. All statistics previously presented for RDA also apply for partial RDA.
-> NOTE : It is really important to note that while using "...", to avoid any confusion, we will have to pass arguments by name.
+Partial RDA is thus a powerful tool when users what to assess the effect of environmental variables on species composition while taking into account the species variation due to other environmental variables with no interest. It can also be used to control for well-known linear effects, isolate the effect of a single explanatory variable or analyse related samples. In the example below, we will assess the effect of water chemistry on fish species abundances partialling out the effect of physiography.
-  * Allow the user to input an indefinite number of arguments. The value of each argument will then have to be checked manually. This is done by transforming "..." as a list and iterating over it. Let's create a sum function that takes an indefinite number of arguments.
+In R, partial RDA is performed in the same way as RDA using rda() with the addition of a condition term:
-<code rsplus>
-sum2 <- function(...){
-  args <- list(...)
-  result <- 0
-  for (i in args)  {
-    result <- result + i
-  }
-  return (result)
-}
-sum2(2, 3)
+<code rsplus | Partial RDA with rda()>
-sum2(2, 4, 5, 7688, 1)
+#Divide the env2 dataframe into two dataframes:
+envtopo <- env[, c(1:3)] # Physiography : explanatory dataset 1
+names(envtopo)
+envchem <- env[, c(4:10)] # Water quality : explanatory dataset 2
+names(envchem)
+#Run the partial RDA
+spechem.physio=rda(spe.hel, envchem, envtopo)
+summary(spechem.physio, display=NULL)
+#or
+spechem.physio2=rda(spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo
+                    + Condition(alt + pen + deb), data=env)
+#Extract the results
+summary(spechem.physio, display=NULL)
+#Calculate the adjusted R2 of the partial RDA
+(R2adj <- RsquareAdj(spechem.physio)$adj.r.squared)
+#Test the significance of the axes in partial RDA
+anova.cca(spechem.physio, step=1000)
+anova.cca(spechem.physio2, step=1000, by="axis")
+#Construct the triplots
+#Scaling 1
+windows(title="Partial RDA scaling 1")
+plot(spechem.physio, scaling=1, main="Triplot partial RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=1),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(spechem.physio, display="species", choices=c(1), scaling=1),
+       scores(spechem.physio, display="species", choices=c(2), scaling=1),
+       col="black",length=0)
+text(scores(spechem.physio, display="species", choices=c(1), scaling=1),
+     scores(spechem.physio, display="species", choices=c(2), scaling=1),
+     labels=rownames(scores(spechem.physio, display="species", scaling=1)),
+     col="black", cex=0.8)
+arrows(0,0,
+      scores(spechem.physio, display="bp", choices=c(1), scaling=1),
+      scores(spechem.physio, display="bp", choices=c(2), scaling=1),
+      col="red")
+text(scores(spechem.physio, display="bp", choices=c(1), scaling=1)+0.05,
+     scores(spechem.physio, display="bp", choices=c(2), scaling=1)+0.05,
+     labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=1)),
+     col="red", cex=1)
+#Scaling 2
+windows(title="Partial RDA scaling 2")
+plot(spechem.physio, scaling=2, main="Triplot partial RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=2),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(spechem.physio, display="species", choices=c(1), scaling=2),
+       scores(spechem.physio, display="species", choices=c(2), scaling=2),
+       col="black",length=0)
+text(scores(spechem.physio, display="species", choices=c(1), scaling=2),
+     scores(spechem.physio, display="species", choices=c(2), scaling=2),
+     labels=rownames(scores(spechem.physio, display="species", scaling=2)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(spechem.physio, display="bp", choices=c(1), scaling=2),
+       scores(spechem.physio, display="bp", choices=c(2), scaling=2),
+       col="red")
+text(scores(spechem.physio, display="bp", choices=c(1), scaling=2)+0.05,
+     scores(spechem.physio, display="bp", choices=c(2), scaling=2)+0.05,
+     labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=2)),
+     col="red", cex=1)
 </code>
-=== Return value ===
-As shown in the previous example, if we want to be able to save the result of our function and be able to use it later, we have to return it at the end using return(). It is important to note that only one return value can be given by a function. If you want to return more than one object, you will have to use composite object such as lists or dataframes.
+This RDA is significant (p<0.001) as well as the two first canonical axis. Water chemistry explained 31.89% of the variance of fish species composition, physiographic covariables explained 41.53% of this variation and the unexplained variation is 26.59%. The adjusted R² of this RDA is 24.13%. The triplot looks like:
-Also, it is important to note that the function ends once it reaches the return() keyword.
-<code rsplus>
-returntest <- function(a, b) {
+{{ :doubs_rdapart2.png |}}
-  return (a) # The function exits here
-  a <- a + b # Not interpreted
-  return (a + b) # Not interpreted
-}
-returntest(2, 3) # R will by default print the return value of your function
-c <- returntest(2, 3) # to save it, don't forget to assign it to another variable
-c
-</code>
------
+**Challenge 2**
-**Exercise 6**
-Using what you learned so far on functions and flow control, create a function **bigsum** that takes two arguments **a** and **b** and :
+Run the partial RDA of the mite environmental variables of the mite species abundances partialling out for the substrate variables (SubsDens, WaterCont and Substrate). Is the model significant? Which are the significant axes? Interpret the obtained triplot.
-  * returns 0 if the sum of a and b is strictly less than 50
-  * returns the sum of a and b otherwise
-++++ Exercise 6 : Answer|
+**Challenge 2**: Solution
-<file rsplus| Exercise6 Answer>
-bigsum <- function(a, b) {
-  result <- a + b
-  if (result < 50) {
-    return(0)
-  } else {
-    return (result)
-  }
-}
-</file>
-++++
------
-==== Accessibility of variables ====
-When working with flow control structures and functions, it is essential to always keep in mind where your variables are and whether they are defined and accessible.
+<hidden>
-Here are some tips to keep in mind
-  * Variables defined inside a function are not accessible outside
+Here is some potential code:
-  * Variables defined outside a function are accessible inside. But it is **NEVER** a good idea to use them inside as your function might not function anymore if the variable doesn't exist
-<code rsplus>
+<code rsplus | Mite partial RDA>
+#Partial RDA
+mite.spe.subs=rda(mite.spe.hel ~ Shrub + Topo
+                  + Condition(SubsDens + WatrCont + Substrate), data=mite.env)
+#Summary
+summary(mite.spe.subs, display=NULL)
+(R2adj <- RsquareAdj(mite.spe.subs)$adj.r.squared)
+#Significant axes
+anova.cca(mite.spe.subs, step=1000)
+anova.cca(mite.spe.subs, step=1000, by="axis")
+#Triplot scaling 1
+windows(title="Partial RDA scaling 1")
+plot(mite.spe.subs, scaling=1, main="Triplot partial RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(mite.spe.subs, display="sites", choices=c(1,2), scaling=1),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(mite.spe.subs, display="species", choices=c(1), scaling=1),
+       scores(mite.spe.subs, display="species", choices=c(2), scaling=1),
+       col="black",length=0)
+text(scores(mite.spe.subs, display="species", choices=c(1), scaling=1),
+     scores(mite.spe.subs, display="species", choices=c(2), scaling=1),
+     labels=rownames(scores(mite.spe.subs, display="species", scaling=1)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(mite.spe.subs, display="bp", choices=c(1), scaling=1),
+       scores(mite.spe.subs, display="bp", choices=c(2), scaling=1),
+       col="red")
+text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=1)+0.05,
+     scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)+0.05,
+     labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)),
+     col="red", cex=1)
+#Triplot scaling 2
+windows(title="Partial RDA scaling 2")
+plot(mite.spe.subs, scaling=2, main="Triplot partial RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(mite.spe.subs, display="sites", choices=c(1,2), scaling=2),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(mite.spe.subs, display="species", choices=c(1), scaling=2),
+       scores(mite.spe.subs, display="species", choices=c(2), scaling=2),
+       col="black",length=0)
+text(scores(mite.spe.subs, display="species", choices=c(1), scaling=2),
+     scores(mite.spe.subs, display="species", choices=c(2), scaling=2),
+     labels=rownames(scores(mite.spe.subs, display="species", scaling=2)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(mite.spe.subs, display="bp", choices=c(1), scaling=2),
+       scores(mite.spe.subs, display="bp", choices=c(2), scaling=2),
+       col="red")
+text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=2)+0.05,
+     scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)+0.05,
+     labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)),
+     col="red", cex=1)
+</code>
-rm(list=ls()) # first let's remove everything to avoid any confusion
+{{ :partial_rda_ch6.png |}}
-var1 <- 3     # var1 is defined outside our function
+This RDA is significant (p<0.001) as well as the first canonical axis. Environmental variables explained 9.81% of the variance of mite species composition, substrate covariables explained 42.84% of this variation and the unexplained variation is 47.35%. The adjusted R² of this RDA is 8.33%.
-vartest <- function() {
+</hidden>
-  a <- 4      # a is defined inside
+=====2.3 Variation partitioning by partial RDA=====
-  print(a)    # print a
-  print(var1) # print var1
-}
-a             # print a. It doesn't work, a can be seen only inside the function
-vartest()     # calling vartest() will print a and var1
-rm(var1)      # remove var1
-vartest()     # calling the function again doesn't work anymore
-</code>
+Variation partitioning is a type of analysis that combines RDA and partial RDA to divide the variation of a response variable among two, three or four explanatory data sets. Variation partitioning are generally represented by Venn diagram in which the percentage of explained variance by each explanatory data set (or combination of data stets) is reported.
-Instead, use arguments!! Inside a function, arguments names will take over other variable names.
+In the case of two datasets (below):
+- Fraction a + b +c is the explained variance by the two datasets calculated using a RDA of y by X + W. - Fraction d is the unexplained variance by the two datasets calculated using the same RDA as above. - Fraction a is the explained variance by the X data set only calculated using a partial of y by X with W as covariables. - Fraction c is the explained variance by the W data set only calculated using a partial of y by W with X as covariables. - Fraction b is calculated by subtraction, i.e. b = [a + b] + [b + c] - [a + b + c].
-<code rsplus>
+{{ :vennd_varpart.png |}}
-var1 <- 3     # var1 is defined outside our function
-vartest <- function(var1) {
-  print(var1) # print var1
-}
-vartest(8)    # Inside our function var1 is now our argument and takes its value
-var1          # var1 still has the same value
-</code>
-Be very careful when creating variables inside a conditionnal statement as the variable can never be created and cause errors
+Venn diagram of partition of the variation of a response variable y among two sets of explanatory variables X and W (from Legendre and Legendre 2012).
-<code rsplus>
+Variation partitioning is thus an indicated analysis when user what to relate the abundance of species in a community to various type of environmental variables, for example abiotic vs biotic variables, large-scale versus small-scale variables, etc. In the next example, we will partition the variation of fish species composition between chemical and physiographic variables.
-a <- 3
-if (a > 5) {
-  b <- 2
-}
-a + b    # Error! b doesn't exist
-</code>
-Usually it is a good practice to define variables outside the conditions and then modify their value to avoid any problem
+In R, variation partitioning is performed using the function varpart(). Venn diagrams can also be drawn using the function plot().
-<code rsplus>
-a <- 3
-b <- 0
-if (a > 5) {
-  b <- 2
-}
-a + b
-</code>
-===== Good practices and how to speed up your code =====
-==== Good practices ====
-Here are some programming tips that can make your life easier, help achieve greater readability and makes sharing and reusing your code a lot less painful. Having a clear to read code will reduce the time you'll spend to understand it so it's never time lost.
+<code rsplus | varpart()>
+?varpart
-=== Keep a clean and nice code ===
+vegandocs("partitioning.pdf")
-One thing that often helps the most when reading programming code, is to have a nicely formatted code, well spaced and well indented code. Some programming stardards exist to help achieve a great consistency but usually it ends up to one's personnal preferences. Here are some tips that can help:
+#Variation partitioning with all explanatory variables
-  * Use spaces between and after your operators
+spe.part.all <- varpart(spe.hel, envchem, envtopo)
-  * Use consistentely the same assignation operator. `<-` is often preferred, `=` is ok but don't switch all the time between the two
+spe.part.all
-  * Use brackets when using flow control statements, even if it's for one line. Each statement inside brackets should be indented by two spaces. The closing brackets would usually be all by themselves on a seperate line, except when preceding an else statement. This helps greatly when trying to identify where we are, especially with a lot of nested loops/conditions.
+windows(title="Variation partitioning - all variables")
-  * Define each variable on its own line
+plot(spe.part.all, digits=2)
-Here is some hard to read code
-<code rsplus>
-a<-4;b=3
-if(a<b){
-if(a==0)print("a zero") } else {
-if(b==0){print("b zero")} else print(b)}
 </code>
-Here is a little easier-to-read version. It takes more space but it is easier to see the flow of the code.
-<code rsplus>
+The output looks like:
-a <- 4
+{{ :varpart_output.png |}}
-b <- 3
+{{ :varpart_output_venn.png |}}
-if(a < b){
-  if(a == 0) {
-    print("a zero")
-  }
-} else {
-  if(b == 0){
-    print("b zero")
-  } else {
-    print(b)
-  }
-}
-</code>
-There are some coding styles that can be found around the internet. Here is an example that provides a clear code : [[https://google-styleguide.googlecode.com/svn/trunk/Rguide.xml|https://google-styleguide.googlecode.com/svn/trunk/Rguide.xml]]
+In this case, the chemical variables explain 24.10% of the fish species composition, the physiographic variables explain 11.20% of the fish species composition and the interaction of these two types of variables explained 23.30% of the fish species composition. Note that the varpart() function also identify the fractions that can be tested for significance using the function anova.cca().
-=== Use functions whenever possible ===
+Users can also perform variation partitioning between data sets that only contain significant environmental variables:
-Now that you know how to create functions, try to use them everytime you can. Everytime you see some portion of code that is repeated more than two times in your script, you should be thinking "Hmmm... would it not be better to write a function instead?". If only a part of the code change, try thinking of ways to insert them as arguments inside a function instead. This would reduce the number of errors done by copying/pasting and the time needed to correct them.
-Let's modify the example from exercise 3 and suppose that all CO2 uptake from Mississipi were overestimated by 20 and Quebec underestimated by 50. We could write this
+<code rsplus | varpart with significant variables>
+#RDA of chemistry variables
-<code rsplus>
+spe.chem <- rda(spe.hel~., data=envchem)
-for (i in 1:length(CO2[,1])) {
-  if(CO2$Type[i] == "Mississippi") {
+#Select significant chemistry variables
-    CO2$conc[i] <- CO2$conc[i] - 20
+R2a.all.chem <- RsquareAdj(spe.chem)$adj.r.squared
-  }
+ordiR2step(rda(spe.hel~1, data=envchem),
-}
+           scope= formula(spe.chem), direction= "forward", R2scope=TRUE, pstep=1000)
-for (i in 1:length(CO2[,1])) {
+names(envchem)
-  if(CO2$Type[i] == "Quebec") {
+(envchem.pars <- envchem[, c( 4, 6, 7 )])
-    CO2$conc[i] <- CO2$conc[i] + 50
-  }
+#RDA with other environmental variables
-}
+spe.topo <- rda(spe.hel~., data=envtopo)
+R2a.all.topo <- RsquareAdj(spe.topo)$adj.r.squared
+ordiR2step(rda(spe.hel~1, data=envtopo),
+           scope= formula(spe.topo), direction= "forward", R2scope=TRUE, pstep=1000)
+names(envtopo)
+envtopo.pars <- envtopo[, c(1,2)]
+#Varpart
+spe.part <- varpart(spe.hel, envchem.pars, envtopo.pars)
+windows(title="Variation partitioning - parsimonious subsets")
+plot(spe.part, digits=2)
+#Tests of significance
+anova.cca(rda(spe.hel, envchem.pars), step=1000) # Test of fractions [a+b]
+anova.cca(rda(spe.hel, envtopo.pars), step=1000) # Test of fractions [b+c]
+env.pars <- cbind(envchem.pars, envtopo.pars)
+anova.cca(rda(spe.hel, env.pars), step=1000) # Test of fractions [a+b+c]
+anova.cca(rda(spe.hel, envchem.pars, envtopo.pars), step=1000) # Test of fraction [a]
+anova.cca(rda(spe.hel, envtopo.pars, envchem.pars), step=1000) # Test of fraction [c]
 </code>
-Or we could do this instead.
-<code rsplus>
+Now, the chemical variables explain 25.30% of the fish species composition, the physiographic variables explain 14.20% of the fish species composition and the interaction of these two types of variables explained 19.60% of the fish species composition. All these fractions are significant (p<0.001).
-recalibrate <- function(CO2, type, bias) {
-  for (i in 1:nrow(CO2)) {
-    if(CO2$Type[i] == type) {
-      CO2$conc[i] <- CO2$conc[i] + bias
-    }
-  }
-  # we have to return our new dataset because the original is not modified
-  return (CO2)
-}
-newCO2 <- recalibrate(CO2, "Mississipi", -20)
-# Note that we recalibrate our newCO2 dataset here because the original CO2 is not modified
-newCO2 <- recalibrate(newCO2, "Quebec", +50)
-</code>
-And now, we realize that what we modified was not the uptake but the concentration... We now have to change all occurences of //CO2\$conc[i]// by //CO2\$uptake[i]//. In the first case, it means we have to change it 4 times, and only 2 times in the second one! (At this point, you might think that this is not much and you can do it with a simple search/replace and you might be right but this is just a simple example! Imagine you had 10 locations instead of 2. A good programmer is a lazy one. Also, admit it, it looks way cooler with a function...)
+**Challenge 3**: Perform variation partitioning of the mite species abundances with a first dataset for the significant substrate variables (SubsDens, WaterCont and Substrate) and a second dataset for the significant other variables (Shrud and Topo). What proportion of the variation are explained by each dataset? What are the significant fractions ?
+**Challenge 3**: Solution
-=== Give meaningful variable and function names ===
+<hidden>
+This is what your code may look like:
-This help to see at first glance what does what. Be extra careful when choosing your argument names when creating a function since it's what users will see. However it is also good to choose short names to avoid having to type them all the time and making typos so a good balance should be reached.
+<code rsplus | varpart with significant variables>
+str(mite.env)
-Here is what our previous example could look like with vague names. It now requires a little more work to understand what this function does.
+(mite.subs=mite.env[,c(1,2,3)]) #First set of variables outlined in challenge
+(mite.other=mite.env[,c(4,5)]) #Second set of variables outlined in challenge
-<code rsplus>
-rc <- function(c, t, b) {
+#RDA for mite.subs
-  for (i in 1:nrow(c)) {
+rda.mite.subs <- rda(mite.spe.hel~., data=mite.subs)
-    if(c$Type[i] == t) {
+R2a.all.subs <- RsquareAdj(rda.mite.subs)$adj.r.squared
-      c$uptake[i] <- c$uptake[i] + b
-    }
+#Forward selection for mite.subs
-  }
+ordiR2step(rda(mite.spe.hel~1, data=mite.subs),
-  return (c)
+           scope= formula(rda.mite.subs), direction= "forward", R2scope=TRUE, pstep=1000)
-}
+names(mite.subs)
+(mite.subs.pars <- mite.subs[, c(2, 3)])
+#RDA for mite.other
+rda.mite.other <- rda(mite.spe.hel~., data=mite.other)
+R2a.all.other <- RsquareAdj(rda.mite.other)$adj.r.squared
+#Forward selection for mite.other
+ordiR2step(rda(mite.spe.hel~1, data=mite.other),
+           scope= formula(rda.mite.other), direction= "forward", R2scope=TRUE, pstep=1000)
+names(mite.other)
+(mite.other.pars <- mite.other[, c(1,2)])
+#Variation partitioning
+(mite.spe.part <- varpart(mite.spe.hel, ~WatrCont+Substrate, ~Shrub+Topo,
+                          data=mite.env))
+windows(title="Variation partitioning - parsimonious subsets")
+plot(mite.spe.part, digits=2)
+# Tests of all testable fractions
+anova.cca(rda(mite.spe.hel~ WatrCont+Substrate, data=mite.env), step=1000) # Test of fractions [a+b]
+anova.cca(rda(mite.spe.hel~Shrub+Topo, data=mite.env), step=1000) # Test of fractions [b+c]
+(env.pars <- cbind(mite.env[,c(2,3,4,5)]))
+anova.cca(rda(mite.spe.hel~ WatrCont+Substrate+Shrub+Topo, data=env.pars), step=1000) # Test of fractions [a+b+c]
+anova.cca(rda(mite.spe.hel~WatrCont+Substrate + Condition(Shrub+Topo), data=env.pars), step=1000) # Test of fraction [a]
+anova.cca(rda(mite.spe.hel~Shrub+Topo+ Condition(WatrCont+Substrate ), data=env.pars), step=1000) # Test of fraction [c]
 </code>
+In this case, substrate variables explain 14.00% of the mite species composition, the other environmental variables explain 9.1% of the mite species composition and the interaction of these two types of vriables explained 16.90% of the mite species composition. All these fractions are significant (p<0.001).
+</hidden>
-=== Comments ===
-Even with meaningful names, it's never a bad thing to add comment to describe everything your code does, be it the purpose of a function, how to use its arguments or a detailed step by step of the function.
-<code rsplus>
+======3. Multivariate regression tree======
-## recalibrates the CO2 dataset by modifying the CO2 uptake concentration
-## by a fixed amount depending on the region of sampling
-# Arguments
-# CO2: the CO2 dataset
-# type: the type that need to be recalibrated. Values: "Mississippi" or "Quebec"
-# bias: the amount to add to the concentration uptake. Use negative values for overestimations
-recalibrate <- function(CO2, type, bias) {
-  for (i in 1:nrow(CO2)) {
-    if(CO2$Type[i] == type) {
-      CO2$uptake[i] <- CO2$uptake[i] + bias
-    }
-  }
-  # we have to return our new dataset because the original is not modified
-  return (CO2)
-}
-</code>
+Multivariate regression tree (MRT) is a constrained clustering technique. Introduced by De’ath (2002), MRTs allow the partitioning of a quantitative response matrix by a matrix of explanatory variables constraining (guiding) on where to divide the data of the response matrix. RDA and MRT are both regression techniques, the former explaining the global structure of relationships through a linear model, the latter better highlighting local structures and interactions among variables by producing a tree model.
-===== 3. Speeding up your code =====
+Advantages of the MRT compared to the RDA:
+    * does not make assumptions about the shape of the relationships between species and environmental variables (quantitative or categorical),
+    * is robust in dealing with missing values
+    * is robust in dealing with collinearity among the explanatory variables
+    * is insensitive to transformations of the explanatory variables, which allows the use of raw values
+    * the outcome, the tree, is easy to interpret, especially to a non-scientist audience.
+The MRT technique splits the data into clusters of samples similar in their species composition based on environmental value thresholds. It involves two procedures running at the same time: 1) the computation of the constrained partitioning of the data, and 2) the calculation of the relative error of the successive partitioning levels by multiple cross-validations. The function mvpart() from the package mvpart computes both the partition and the cross-validation.
-Here are some programming tips to program more efficiently with R and achieve greater performance and faster code. Note that before optimizing your code, you should always make sure to have a working code beforehand. A slow working code is always better than a fast broken implementation. Also, sometimes, optimizing is just not the way to go. If you spend 2 hours rewritting code to gain a few seconds, then it might just not be worth it...
+A quick note on MRT terminology:
-==== Before we start : profiling our code ====
+Leaf: Terminal group of sites
-If we want to optimize our code, we will need to know how much time each task takes to perform.
+Node: Point where the data splits into two groups. It is characterized by a threshold value of an explanatory variable.
-The simplest way to do that is to use the function system.time(//expression//)
+Branch: Each group formed by a split
-<code rsplus>
-system.time({
-a <- 0
-  for (i in 1:1000) {
-    a <-  a + i
-  }
-})
-</code>
-Note that most of the time R works really fast and you will need to have some heavy computing to do or your time might not even be recorded. It is usually recommended to perform iterations of the task you want to profile or work with really big datasets.
+**1- Constrained partitioning of the data**
-<code rsplus>
+First, the method computes all possible partitions of the sites into two groups. For each quantitative explanatory variable, the sites will be sorted in the ascending values of the variables; for categorical variables, the sites will be aggregated by levels to test all combinations of levels. The method will split the data after the first object, the second object and so on, and compute the sum of within-group sum of squared distances to the group mean (within-group SS) for the response data. The method will retain the partition into two groups minimizing the within-group SS and the threshold value/level of the explanatory variable. These steps will be repeated within the two subgroups formed previously, until all objects form their own group. In other words, when each leaf of the tree contains one object.
-system.time(replicate(1000, {
-  a <- 0
-  for (i in 1:1000) {
-    a <-  a + i
-  }
-}))
-</code>
-One other simple and useful tool is the function Rprof(). The main advantage of Rprof() is that it saves information about time spent in each function in a file that you can access later. Here's how to use it.
+**2- Cross-validation and pruning the tree**
-<code rsplus>
+The mvpart function also performs a cross-validation and identifies the best predictive tree. The cross-validation procedure consists in using a subset of the objects to construct the tree, and to allocate the remaining objects to the groups. In a good predictive tree, objects are assigned to the appropriate groups. The cross-validated relative error (CVRE) is the measure of the predictive error. Without cross-validation, one would retain the number of partitions minimizing the variance not explained by the tree (i.e. the relative error: the sum of the within-group SS over all leaves divided by the overall SS of the data). This is the solution maximizing the R2 so to speak. This approach is explanatory rather than predictive.
-Rprof("profile.txt")  # you can change profile.txt by the filename you want
-for (i in 1:1000) {
-    a <- 0
-    for (i in 1:1000) {
-      a <-  a + i
-    }
-  }
-Rprof()               # This ends the profiling
-summaryRprof("profile.txt")  # Use the filename previously recorded to display the summary of the tasks
-</code>
-Finally, if you want to compare the efficiency or several functions, a very good tool is the package microbenchmark
+Let’s create a multivariate regression tree on the Doubs data.
-<code rsplus>
-install.packages("microbenchmark")
-library(microbenchmark)
-f1 <- function() {
+<code rsplus | mvpart()>
-  a <- 0
+?mvpart
-  for (i in 1:1000) {
-    a <-  a + i
+#Prepare the data: remove “distance from source”
-  }
+env <- subset(env, select = -das)
-}
-microbenchmark(f1(), times=1000) # the argument times allow us to determine how many iterations we want
+# Create the regression tree
+doubs.mrt <- mvpart(as.matrix(spe.hel) ~. ,env,
+            legend=FALSE, margin=0.01, cp=0, xv="pick",
+            xval=nrow(spe.hel), xvmult=100, which=4)
 </code>
-==== First step : thinking a bit! ====
-When you look at your code, often you will realize that there are simpler, way more efficient ways to do what you want and that some operations can be easily removed for added speed.
+At this point, you will need to select the tree who's size (number of groups) is appropriate to the aim of your study from the following graph. This step requires the argument xv="pick". In other words, you must prune the tree by picking the best-fit tree. Indeed, a fully resolved tree is not the desirable outcome. Instead, one is usually interested in a tree including only informative partitions/groups. It is possible to have an a-priori idea of the number of potential groups to be retained as well.
-For example, let's create a function that takes a number **a**. We will add **a** to every number from 1 to 100, and if **a** is less than 5, then we will add 2*a instead. Then we will calculate the sum of all the elements of the sequence.
-Here's a way to do it
+{{ :cross_validation.png |}}
-<code rsplus>
+The graph shows the relative error RE (in green) and the cross-validated relative error CVRE (in blue) of trees of increasing size. The red dot indicates the solution with the smallest CVRE, and the orange dot shows the smallest tree within one standard error of CVRE. It has been suggested that instead of choosing the solution minimizing CVRE, it would be more parsimonious to opt for the smallest tree for which the CVRE is within one standard error of the tree with the lowest CVRE (Breiman et al. 1984). The green bars at the top indicate the number of times each size was chosen during the cross-validation process.
-f2 <- function(a) {
-  # initialize our result
-  result <- 0
-  # iterate on the sequence from 1 to 100
-  for (i in 1:100) {
-    if (a < 5) {
-      # a is < 5, we add 2*a to the sequence element. We save it in result
-      result <- result + i + (2*a)
-    } else {
-      # a is >= 5, we add only a
-      result <- result + i + a
-    }
-  }
-  return(result)
-}
-f2(4)
-</code>
-This is ok and does what we want. However we have a lot of useless steps in our code. For example, we don't really need to perform our condition on each iteration, as the result will always be the same.
+This graph is interactive, which means you will have to click on the blue point corresponding your choice of tree size. Once you do so, the corresponding multivariate regression tree will appear. If you click on the orange dot, the following tree appears.
-So let's take it out of the loop
-<code rsplus>
-f3 <- function(a) {
-  # initialize our result
-  result <- 0
-  # Check if a < 5. If true, a becomes 2*a
-  if (a < 5) {
-   a <- 2 * a
-  }
-  # We don't even need an else here since a remains the same otherwise
-  # iterate on the sequence from 1 to n
-  for (i in 1:100) {
-      result <- result + i + a
-  }
-  return(result)
-}
-f3(4)
+{{ :mrt_1se.png |}}
-microbenchmark(f2(4),
-               f3(4), times=1000)
-</code>
+The statistics at the bottom of the figure are: the residual error (the reciprocal of the R2 of the model, in this case 43.7%), the cross-validated error, and the standard error. This tree has only two leaves separated by one node. This node splits the data into two groups at the threshold altitude value of 361.5m.
-This is just a simple modification but here we sped up our code by almost 40% (results may vary depending on computers). Moreover, our code is easier to read and to understand. Sometimes, we can gain speed and easier code just by thinking a little more of where our conditional statements can be and what they test.
+Each leaf is characterized by a small barplot showing the abundances of the species, its number of sites and its relative error.
-But, using the strengths of R, we can do even better!
+We can compare this tree with the 10-group solution, as suggested by the CVRE criterion, or choose a solution in between, e.g. with 4 leaves to compare.
-<code rsplus>
-f4 <- function(a, n) {
-  result <- 0
-  if (a < 5) {
-    a <- a + 1
-  }
-  result <- sum(1:n + a)
-  return(result)
-}
-f4(4)
+<code rsplus | Compare trees>
-microbenchmark(f3(4),
+# Using the CVRE criterion
-               f4(4), times=1000)
+doubs.mrt.cvre <- mvpart(as.matrix(spe.hel)~., env,
+                 legend=FALSE, margin=0.01, cp=0,xv="pick",
+                 xval=nrow(spe.hel), xvmult=100,which=4)
+# Choosing ourself the best number of partitions
+doubs.mrt.4 <- mvpart(as.matrix(spe.hel)~., env,
+              legend=FALSE, margin=0.01, cp=0, xv="pick",
+              xval=nrow(spe.hel), xvmult=100,which=4)
 </code>
-Wow, here the improvement is much more efficient... But how exactly did this happen? This leads us to our next point
+{{ :mrt_cvre.png |}}
+{{ :mrt_4.png |}}
-==== Vectorization ====
+The 10-group solution has a high EXPLANATORY power but its predictive power (indicated by the cross-validated error) is just slightly better than that of the 2-group solution. The 4-group solution seems to be a good compromise.
-This part is a reminder of what you probably already learnt in the first workshops. However, this is usually overlooked and forgotten and the bad performance of R code can often be attributed to bad vectorization.
-R is meant to function with vectors and many functions in R are optimized for vectorization. To understand that, it is first important to understand how R works at a low level. R is an interpreted language, which means that when you execute R code, you are in reality sending your code to functions programmed in another language (the C language). This slows down the execution of your programs since first you have to interpret the R code then pass it on to other functions. When you create a loop in R, you have to decode every iteration and pass it on.
-Vectorized functions on the other hand are functions that perform operations directly on a vector. Essentially, they also perform a loop on your vector, but the big difference is that they do it in C, not R and are therefore much faster. sum() is one of these functions.
-One of the biggest challenge with R is learning to think and program with vectors and not with single elements. For example, most of the basic operations can be done on vectors.
-<code rsplus>
+More information can be obtained by looking at the summary output.
-v1 <- 1:5
-v2 <- 2:6
-v3 <- 1:3
-v1 + 2      # Addition on a vector : adds 2 to all elements
-v1 + v2     # Adds each element of v2 to v
-v1 + v3     # v1 and v3 are not the same length, then we add from the start of v3 again
-sum(v1)     # Adds all elements of v1 together
-sum(v1, v2) # Sums all elements of v1 and v2
-mean(v1)    # Average of elements in v1
-mean(c(v1, v2)) # Average of elements of v1 and v2. Unlike sum, we have to combine them beforehand
-</code>
-=== Subsetting ===
-To vectorize efficiently, it is also important to be able to extract values from our data quickly.
+<code rsplus | MRT summary>
-To apply treatment on specific elements of a vector or a dataframe, R offers a subsetting tool that can be sometimes way more efficient and easy to write than a mix of loops and conditions.
+summary(doubs.mrt)
-Subsetting is done via the **[** and the **$** (for a dataframe) operators. We can insert directly our conditions inside the **[]** part of the subsetting to quickly extract values from our data. We can also use the function **which()** to test for a condition. which() returns the indexes of the elements that matches the condition.
-<code rsplus>
-v1 <- 1:10
-v1[7]      # Extracts the 7th value
-v1[v1 > 5] # Extracts values > 5 only
-v1[which(v1 > 5)]  # same as before
 </code>
-In dataframes, **$** allows to acces a column by name. We can also do this by providing the name of the column directly
+{{ :doubs_mrt_summary.png |}}
-<code rsplus>
+CP stands for “complexity parameter”, which is the equivalent of the variance explained by each node. The CP at nsplit 0 is the R2 of the whole tree. The summary then outlines, for each node, the best threshold values to split the data. While informative, this output is very dense. A more detailed and yet more manageable output can be generated by using the wrapper from the function MRT() of the MVPARTwrap package.
-CO2 <- read.csv("co2_good.csv")
+Plus, this other function allows identification of discriminant species.
-CO2$Type  # Prints columns Type
-CO2[, "Type"] # Same as above
-CO2[CO2$Type == "Quebec", ] #Extracts all rows of the CO2 dataset where the Type is "Quebec"
-</code>
------
+<code rsplus | Find discriminant and indicator species>
-**Exercise 7**
+# Find discriminant species with MRT results
+doubs.mrt.wrap<-MRT(doubs.mrt,percent=10,species=colnames(spe.hel))
+summary(doubs.mrt.wrap)
-Create a new function recalibrate2() that is a rewrite the function recalibrate seen earlier using subsetting and vectorization techniques. The new function should not be longer than 3 lines.\\
+# Extract indval p-values
-Reminder:
+doubs.mrt.indval<-indval(spe.hel,doubs.mrt$where)
+doubs.mrt.indval$pval
-<code rsplus>
+# Extract indicator species of each node, with its indval
-recalibrate <- function(CO2, type, bias) {
+doubs.mrt.indval$maxcls[which(doubs.mrt.indval$pval<=0.05)]
-  for (i in 1:nrow(CO2)) {
+doubs.mrt.indval$indcls[which(doubs.mrt.indval$pval<=0.05)]
-    if(CO2$Type[i] == type) {
-      CO2$uptake[i] <- CO2$uptake[i] + bias
-    }
-  }
-  return (CO2)
-}
 </code>
-++++ Exercise 7 : Answer|
+{{ ::doubs_mrt_discriminant.png |}}
-<code rsplus>
+{{ ::doubs_mrt_finalpart.png |}}
-recalibrate2 <- function(CO2, type, bias) {
-  # First get the indexes of the data with the good type
-  # Thinking tip : since we use the indexes twice below, instead of using which() twice, let's
-  # do it only once and save the result!
-  idx <- which(CO2$Type == type)
-  # Modify only the data concerned using indexes.
-  CO2$uptake[idx] <- CO2$uptake[idx] + bias
-  return (CO2)
-}
-# Check the results are the same
+The main discriminant species of the first split are TRU, VAI and ABL. TRU and VAI contribute highly to the left leaf, and ABL is the most indicative species of the sites at lower altitude (<361.5m). This output also indicates which sites are included in each leaf.
-all.equal(recalibrate(CO2, "Quebec", 20), recalibrate2(CO2, "Quebec", 20))
-# Check that this is indeed way faster
-microbenchmark(recalibrate(CO2, "Quebec", 20),
-               recalibrate2(CO2, "Quebec", 20))
-</code>
-++++
------
-==== Growing objects ====
+{{ ::doubs_mrt_indval.png |}}
-Vectorizing is good, but sometimes it is hard to do and it might take you more time to do than just do a simple loop. Sometimes you just need loops and you should use them when you can. But when using them, one thing you should be aware of if you want decent performance are growing objects. That is objects that are getting bigger and bigger with each iteration.
+The second part of the code allows us to test the significance of the indicator value of each species through a permutation test. For each significant indicator species, we extracted the leaf number and the indicator value. In this particular case, TRU, VAI and LOC are all significant species of the left leaf, TRU having the highest indicator value (0.867).
-Let's illustrate this really simply by creating a function that iterates over a sequence and create a vector with it. We will compare 2 ways of doing that : by growing our object and by preallocating our result and just modify it.
-<code rsplus>
+**Challenge 4**: Run the multivariate regression tree for the mite data. Select the minimum size of tree within one SE of the CVRE. What is the proportion of variance explained by this tree? How many leaves contain this tree? What are the discriminant species?
-growing <- function(n) {
-  # declare our result
-  result <- NULL
-  for (i in 1:n) {
-    # create our result by growing our object
-    result <- c(result, i)
-  }
-  return(result)
-}
-growing2 <- function(n) {
+**Challenge 4** - Solution
-  # declare our result : here we create a vector of length n with 0 in it
-  result <- numeric(n)
-  for (i in 1:n) {
-    # now we just modify our value instead of recreating the vector
-    result[i] <- i
-  }
-  return(result)
-}
-</code>
-Now let's compare their speeds
+<hidden>
-<code rsplus>
+<code rsplus | Create a MRT with the mite data>
-system.time({
+mite.mrt<-mvpart(data.matrix(mite.spe.hel)~.,mite.env,
-  growing(10000)
+legend=FALSE,margin=0.01,cp=0,xv="pick",
-})
+xval=nrow(mite.spe.hel),xvmult=100,which=4)
-system.time({
+summary(mite.mrt)
-  growing2(10000)
-})
-</code>
-With a vector of 10000 elements, speeds are still comparable and under the second. Now let's just use 50000 elements
-<code rsplus>
-system.time({
-  growing(50000)
-})
-system.time({
-  growing2(50000)
-})
-</code>
-By multiplying by only 5, it suddenly takes us several seconds just to create a single vector while modifying  preallocated vector is still almost instantaneous. What happened here? The reason is that when you call a function, arguments are first copied before being passed on to your function. So when you write // result <- c(result, i) // each time, **result** is copied before being passed on to c(). As result grows which each iteration, in each iteration it takes more and more time to copy it. The bigger the end object, the more time it will take. This is why it is always better to create your result object before your loop if you already know what size it will be.
-This is especially valid when working with dataframes and functions such as **rbind()** and **cbind()**.
-Unfortunately, preallocating dataframes does not work that well and there is a better way, albeit a little more complicated. It is to store each row (or column) in a preallocated list first and then call **rbind()** (or **cbind()**)on all elements via the function **do.call()**. The function **do.call()** allows you to execute a named function on a list of arguments. This way **rbind()** is called only once, at the end, which removes the problem of copying the growing object each time we call it.
-<code rsplus>
-growingdf <- function(n, row) {
+mite.mrt.wrap<-MRT(mite.mrt,percent=10,species=colnames(mite.spe.hel))
-  # preallocate our dataframe
+summary(mite.mrt.wrap)
-  df <- data.frame(numeric(n), character(n), stringsAsFactors=FALSE)
-  for (i in 1:n) {
-    # replace the ith row with row
-    df[i,] <- row
-  }
-  return(df)
-}
-growingdf2 <- function(n, row) {
+mite.mrt.indval<-indval(mite.spe.hel,mite.mrt$where)
-  # this is the way to allocate a list with n elements
+mite.mrt.indval$pval
-  df <- vector("list", n)
-  for (i in 1:n) {
-    # put row in the ith element
-    df[[i]] <- row
-  }
-  return(do.call(rbind, df))
-}
-# store our row in a list since we have different types
+mite.mrt.indval$maxcls[which(mite.mrt.indval$pval<=0.05)]
-row <- list(1, "Hello World")
+mite.mrt.indval$indcls[which(mite.mrt.indval$pval<=0.05)]
-microbenchmark(growingdf(5000, row),
-               growingdf2(5000, row),
-               times=10)
 </code>
+.6% of the variation in the mite species assemblage across sites is explained by the partition of the sites based on water content of the substrate (at 385.1 mg/l). LCIL is a discriminant species of sites with higher water content, and has an indicator value of 0.715.
-==== The apply family ====
+</hidden>
-To prevent the problem of growing objects in loops and to facilitate the application of functions on objects like dataframes, R offers us what we will call the apply functions (because they all have apply in their name...). It is a group of functions that will execute another function on a given object type. Their use only differ depending on the type of object the function is applied to or the type of the return value.
-The apply functions are not always the best choice performance wise as they will usually hide a for loop written in R in their code. However, they can greatly reduce the programming time needed by the ease-of-use they provide.
+======4. Linear discriminant analysis======
-One of the most popular is simply apply() that executes a function on the rows or columns of a dataframe or a matrix. This function takes 3 main arguments:
+Linear discriminant analysis (LDA) is a constrained (canonical) technique that allows you to determine how well your independent set of variables explains an a priori grouping. This grouping may have been obtained from a previous clustering analysis (see Workshop 8) or from a hypothesis (e.g. grouping is based on sites at different latitudes or different treatments). An LDA can also be used to classify new data into these pre-determined groups. You can imagine some useful applications of this technique including assessing which population a fish should be classified in based on morphology or classifying whether a new paper is a freshwater, marine or terrestrial study based on the abstract of papers in those pre-determined biomes.
-  * the objet on which we want to apply the function
-  * the margin or the subscript on which we want to apply the function. 1 is for rows, 2 is for columns
-  * the function to apply
-  * the eventual arguments to the function supplied
-<code rsplus>
+LDA computes discriminant functions from standardized descriptors. These coefficients quantify the relative contributions of the (standardized) explanatory variables to the discrimination of objects. Identification functions can be computed from the original (not standardized) descriptors to classify new data into pre-determined groups.
-df <- data.frame(1:100, 101:200)
+Let’s continue to work with the Doubs fish data. First we must ensure that the within-group covariance matrices of the explanatory variables are homogeneous – a condition necessary for the application of LDA.
-# Sum on rows
-apply(df, 1, sum)
-# Mean on columns
-apply(df, 2, mean)
-# we can also supply additionnal arguments to the function
-apply(df, 2, mean, na.rm=TRUE)
-# we can also define a function directly. The first argument is always what
-# we iterate on. Here each row is treated as a vector of numbers, as we can see with
-# the str() function
-apply(df, 1, function(x){str(x)})
-# We can also add other arguments
-apply(df, 1, function(x, y){x[2] - x[1] + y}, y=5)
-</code>
-All apply functions work on the same model. Performance wise, the most interesting are probably lapply and vapply since they are primitive function written in C. lapply returns a list of the same length as the original object. vapply allows you to specify the format of the value returned by your function. It can be a vector or an array.
+First we want to make an a priori classification that is independent from the environmental data set. We know that there is a general relationship that indicates environmental variables change with latitude (Budyko 1969). Here we will classify our Doubs fish sites based on latitude to determine how well the environmental factors explain our latitude grouping. Our groups are determined by simply dividing the range of latitudes equally into three groups and then assigning each site to a group depending on where they fall along the divided range.
-<code rsplus>
+<code rsplus | Load spatial data and classify Doubs fish sites based on latitude>
-a <- list(1:100, 101:200)
+#load spatial data to determine groups
-# apply mean to each element of the list
+spa <- read.csv ('http://www.davidzeleny.net/anadat-r/data-download/DoubsSpa.csv', row.names = 1)
-lapply(a, mean)  # we get a list as a result
+spa <- spa[,-8]
-unlist(lapply(a, mean)) # use unlist to get a vector instead
-vapply(a, mean, 0) # the result of mean is a single number, we tell vapply our result will be a number
-</code>
-======= 4. Quick introduction to useful packages in R =======
+#View spatial data
+View (spa)
+#add site numbers
+numbers<-(1:30)
+numbers<-numbers[!numbers%in%8]
+spa$site<-numbers
-===== Knitr =====
+#make groups based on lattitude y<82=group1, 82<y<156=group2, y>156=group3
-Knitr is a package that can be used to generate dynamic reports or web pages from R code. The code is evaluated at the moment the report is generated.
+spa.group<-ddply(.data=spa, .variables=.(x, y, site), .fun= summarise, group = if(y <= 82) 1 else if (y <= 156) 2 else 3)
-Code can be easily written in RStudio use the Markdown language. :
+#order by site
+spa.group<-spa.group[with(spa.group, order(site)), ]
+</code>
-++++ Example Markdown code |
-<code rsplus>
----
-title: "Interesting packages for R / QCBS Workshop on R programming"
-output: html_document
----
-### Data table
+Generally, we would first want to check that the within-group covariance matrices of the explanatory variables are homogeneous by verifying multivariate homogeneity of within-group covariance (MHV). For the purposes of this workshop we will by pass it but more information can be found in Borcard et al. (2011).
-A package to facilitate and to improve the efficiency of certain operations in R.
-```{r}
-library(data.table)
-mydf=data.frame(a=rep(LETTERS,each=1e5),b=rnorm(26*1e5))
-mydt=data.table(mydf)
-setkey(mydt,a) # We set the column that will be used as a key for the data table
-```
+Once we run the LDA we can use the result object to determine 1. What groups the sites are classified in based on the environmental data. 2. What are the posterior probabilities of that the sites to belong to the groups. 3. The percentage of correct classification based on our latitudinal grouping.
-Returns all rows with column a (the key) equal to F
-````{r}
-mydt['F']
-````
-Gives the mean value of column b for each letter in column a.
-```{r}
-mydt[,mean(b),by=a]
-```
-Compare
+<code rsplus | Run the linear discriminant analysis (LDA)>
-```{r}
+#run LDA
-system.time(t1<-mydt[,mean(b),by=a])
+LDA<-lda(env,spa.group[,4])
-```
-### With tapply()
+#classification of the objects based on LDA
-```{r}
+spe.class <- predict(LDA)$class
-system.time(t2<-tapply(mydf$b,mydf$a,mean))
-```
-### With reshape2
+#posterior probabilities of the objects to belong to the groups
-```{r message=FALSE}
+spe.post <- predict(LDA)$posterior
-library(reshape2)
-meltdf=melt(mydf)
-system.time(t3<-dcast(meltdf,a~variable,mean))
-```
-### With plyr
+#table of prior versus predicted classifications
-```{r}
+spe.table <- table(spa.group[,4], spe.class)
-library(plyr)
-system.time(t4<-ddply(mydf,.(a),summarize,mean(b)))
-```
-### With dplyr
+#proportion of correct classification
-```{r message=FALSE}
+diag(prop.table(spe.table, 1))
-library(dplyr)
+</code>
-ti1<-proc.time()
-groups <- group_by(mydf, a)
-t4b <- summarise(groups, total = mean(b))
-eltime<-proc.time()-ti1
-eltime
-```
-### With sqldf
+{{ :lda_spetable.png?300 |}}
-```{r}
-library(sqldf)
-system.time(t5<-sqldf('SELECT a, avg(b) FROM mydf GROUP BY a'))
-```
-### With a for loop
+The results suggest that the environmental factors explain the first, lower latitude, group and the group 3 perfectly but only 83% of the group 2 sites were predicted correctly. What does that tell us about our classification? Perhaps there are stronger delineations in the lower and higher latitude and the group 2 is a mix of both?
-```{r}
-ti1<-proc.time()
-t6<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-for (i in t6$letter ){
-  t6[t6$letter==i,2]=mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-```
+Now what we have some new sites and we want to classify them based on the relationship we have established between our latitudinal grouping and environmental factors using the LDA. Using the predict() function we can load in a new matrix with sites and classify them using the LDA object.
-### With a parallelized FOR loop
+Load in the classifyme.csv file, which contains dummy data from 5 new sites.
-```{r}
-library(foreach)
-library(doMC)
-registerDoMC(4) #Four-core processor
-ti1<-proc.time()
-t7<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-t7[,2] <- foreach(i=t7$letter, .combine='c') %dopar% {
- mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-```
+<code rsplus | Load the classify me data and predict grouping of new data>
+#predicting classification of new data
+#read in new sites
+classify.me<-read.csv("classifyme.csv", header = T)
-### RgoogleMaps
+#predict grouping of new data
-```{r message=FALSE}
+predict.group<-predict(LDA, newdata=classify.me)
-library(RgoogleMaps)
-myhome=getGeoCode('Olympic Stadium, Montreal');
-mymap<-GetMap(center=myhome, zoom=14)
-PlotOnStaticMap(mymap,lat=myhome['lat'],lon=myhome['lon'],cex=5,pch=10,lwd=3,col=c('red'));
-```
+#give classification for each new site
+group.new<-predict.group$class
+</code>
-### Taxize
+{{ :lda_newgroups.png?200 |}}
-```{r message=FALSE, warning=FALSE}
-library(taxize)
-spp<-tax_name(query=c("american beaver"),get="species", db = 'ncbi')
-spp
-fam<-tax_name(query=c("american beaver"),get="family", db = 'ncbi')
-fam
-correctname <- tnrs(c("fraxinus americanus"))
-correctname
-cla<-classification("acer rubrum", db = 'itis')
-cla
-```
-### spocc
+Our new sites, in order, have been classified in groups 1,1, 1, 3 and 3 respectively.
-```{r message=FALSE}
-library(spocc)
-occ_data <- occ(query = 'Acer nigrum', from = 'gbif')
-mapggplot(occ_data)
-```
+**Challenge 5**: Run an LDA for the mite env data (only first two vars) based on four latitudinal groups you create from the mite.xy data set. What group was group 2 most incorrectly grouped into? What proportion of sites was correctly classified in group 1? group 2?
+**Challenge 5**: Solution
-### Combine spocc and RgoogleMaps
+<hidden>
-```{r message=FALSE, warning=FALSE}
+<code rsplus | LDA on mite data>
-occ_data <- occ(query = 'Puma concolor', from = 'gbif')
+mite.xy$site<-seq(1:70)
-occ_data_df=occ2df(occ_data)
+(max(mite.xy[,2])-min(mite.xy[,2]))/4
-occ_data_df<-subset(occ_data_df,!is.na(latitude) & latitude!=0)
-mymap<-GetMap(center=c(mean(occ_data_df$latitude),mean(occ_data_df$longitude)), zoom=2)
-PlotOnStaticMap(mymap,lat=occ_data_df$latitude,lon=occ_data_df$longitude,cex=1,pch=16,lwd=3,col=c('red'));
-```
+mite.xy.group<-ddply(.data=mite.xy, .variables=.(x, y, site), .fun= summarise, group = if(y <= 2.5) 1 else if (y <= 4.9) 2 else if (y <= 7.3) 3 else 4)
+mite.xy.group<-mite.xy.group[with(mite.xy.group, order(site)), ]
-### geonames
+LDA.mite<-lda(mite.env[,1:2],mite.xy.group[,4])
-```{r message=FALSE}
+mite.class <- predict(LDA.mite)$class
-library(geonames)
+mite.post <- predict(LDA.mite)$posterior
-options(geonamesUsername="glaroc")
+mite.table <- table(mite.xy.group[,4], mite.class)
-res<-GNsearch(q="Mont Saint-Hilaire")
+diag(prop.table(mite.table, 1))
-res[,c('toponymName','fclName')]
-dc<-GNcities(45.4, -73.55, 45.7, -73.6, lang = "en", maxRows = 10)
-dc[,c('toponymName')]
-```
 </code>
-++++
+</hidden>
-[[http://qcbs.ca/wp-content/uploads/2014/12/QCBS_Workshop8_knitr.html|View the resulting web page]].
-===== Data Table =====
-[[http://cran.r-project.org/web/packages/data.table/index.html|Data table]] is a very useful package which can facilitate and improve the efficiency of certain operations in R. Data tables are just like data frames. You can even create them from data frames.
-[[http://cran.r-project.org/web/packages/data.table/vignettes/datatable-intro.pdf|Introduction to Data table (PDF)]]
+======5. Some other useful ordination methods======
-<code rsplus>
+<code rsplus | other methods>
-install.packages('data.table')
+?cca #(constrained correspondence analysis)
-library(data.table)
+# Constrained Correspondence Analysis (CCA) is a canonical ordination method similar to RDA that preserve
+# Chi-square distances among object (instead of Euclidean distances in RDA). This method is well suited for the
+# analysis of large ecological gradients.
+?CCorA # Canonical Correlation Analysis
+# Canonical Correlation Analysis (CCorA) differs from RDA given that the two matrices are considered symmetric
+# while in RDA the Y matrix is dependent on the X matrix. The main use of this technique is to test the
+# significance of the correlation between two multidimensional data sets, then explore the structure of the data by
+# computing the correlations (which are the square roots of the CCorA eigenvalues) that can be found between
+# linear functions of two groups of descriptors.
+help(coinertia, package=ade4) # Coinertia Analysis
+#Coinertia Analysis (CoIA) is a symmetric canonical ordination method that is appropriate to compare pairs
+# of data sets that play equivalent roles in the analysis. The method finds a common space onto which the objects
+# and variables of these data sets can be projected and compared. Compared to CCorA, co-inertia analysis
+# imposes no constraint regarding the number of variables in the two sets, so that it can be used to compare
+# ecological communities even when they are species-rich. Co-inertia analysis is not well-suited, however, to
+# analyse pairs of data sets that contain the same variables, because the analysis does not establish one-to-one
+# correspondences between variables in the two data sets; the method does not ‘know’ that the first variable is the
+# same in the first and the second data sets, and likewise for the other variables.
+help(mfa, package=ade4) # Multiple Factorial Analysis
+# Multiple factor analysis (MFA) can be used to compare several data sets describing the same objects. MFA
+# consists in projecting objects and variables of two or more data sets on a global PCA, computed from all data
+# sets, in which the sets receive equal weights.
+# Spatial analysis can be performed using packages AEM and PCNM : http://r-forge.r-project.org/R/?group_id=195
 </code>
-Generate very long data frame with one column with letters, and one column with random numbers
-<file rsplus>
-mydf<-data.frame(a=rep(LETTERS,each=1e5),b=rnorm(26*1e5))
-</file>
-Convert the data frame to a data table format.
+======References======
-<file rsplus>
-mydt<-data.table(mydf)
-</file>
-Each data table has to be assigned a key, which is one (or more) of the columns from the table. This key defines the basis for the organization and the sorting of the table.
+Alday & Marrs (2014). A simple test for alternative states in ecological restoration: the use of principal response curves. Journal of Vegetation Science, 17, 302-311.
-<file rsplus>
-setkey(mydt,a)
-</file>
-Once the key is set, we can return all rows with column a (the key) equal to F
+Borcard, Gillet & Legendre (2011). Numerical Ecology with R. Springer New York.
-<file rsplus>
-mydt['F']
-</file>
-Gives the mean value of column b for each letter in column a.
+Breiman, L., J. H. Friedman, et al. (1984). Classification and Regression Trees. Belmont, California, USA, Wadsworth International Group.
-<file rsplus>
-mydt[,mean(b),by=a]
-</file>
-Let's compare the performance of Data table with other methods to achieve the same thing.
+Budyko, M.I. (1969) The effect of solar radiation variations on the climate of the Earth. Tellus, 21(5), 611-619.
-<file rsplus>
-system.time(t1<-mydt[,mean(b),by=a])
-</file>
-**With tapply()**
+Clarke & Warwick (2001). Change in Marine Communities: An Approach to Statistical Analysis and Interpretation 2nd edition. Primer-E Ltd.
-<file rsplus>
-system.time(t2<-tapply(mydf$b,mydf$a,mean))
-</file>
-**With [[http://cran.r-project.org/web/packages/reshape2/index.html|reshape2]]**
-**NOTE**: plyr and reshape2 where covered in [[r_workshop4|Workshop 4]].
-<file rsplus>
-library(reshape2)
-meltdf<-melt(mydf)
-system.time(t3<-dcast(meltdf,a~variable,mean))
-</file>
-**With [[http://cran.r-project.org/web/packages/plyr/index.html|plyr]]**
-, a set of tools to split up a data into homogeneous pieces, apply a function to each piece and combine all the results back together.
-<file rsplus>
-library(plyr)
-system.time(t4<-ddply(mydf,.(a),summarize,mean(b)))
-</file>
-**With [[http://cran.rstudio.com/web/packages/dplyr/vignettes/introduction.html|dplyr]]**
+De'ath, G. (2002). Multivariate regression trees : a new technique for modeling species-environment relationships. Ecology, 83(4), 1105–1117.
-, a newer and faster version of plyr, which is adapted to work only on data frames.
-<file rsplus>
-library(dplyr)
-ti1<-proc.time()
-groups <- group_by(mydf, a)
-t5 <- summarise(groups, total = mean(b))
-eltime<-proc.time()-ti1
-</file>
-**With [[http://cran.r-project.org/web/packages/sqldf/index.html|sqldf]]**. This package allows one to write Structured Query Language commands to perform queries on data frames.
+Gotelli & Ellison (2004). A Primer of Ecological Statistics. Sinaeuer Associates Inc., Sunderland MA.
-<file rsplus>
-library(sqldf)
-system.time(t6<-sqldf('SELECT a, avg(b) FROM mydf GROUP BY a'))
-</file>
-**With a basic FOR loop**
+Legendre & Legendre (2012). Numerical Ecology 3rd edition. Elsevier Science BV, Amsterdam.
-<file rsplus>
-ti1<-proc.time()
-# Initialize an empty data frame with two columns and 26 rows
-t7<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-for (i in t6$letter ){
-  t7[t7$letter==i,2]=mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-</file>
-**With a parallelized FOR loop**
+Poulin, Andersen & Rochefort (2013) A new approach for tracking vegetation change after restoration: a case study with peatlands. Restoration Ecology, 21, 363-371.
-Use [[http://cran.r-project.org/web/packages/foreach/index.html|foreach]] and [[http://cran.r-project.org/web/packages/doMC/index.html|doMC]] packages to run sections of code in parallel on computers with multiple cores. This is particularly suited to speed up some calculations involving FOR loops in which every iteration can be run independently of other iterations. Note that the doMC package may not work on computers running Windows. It should work on Linux or Mac OSX.
-<file rsplus>
-library(foreach)
-library(doMC)
-registerDoMC(4) #Four-core processor
-ti1<-proc.time()
-t8<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-t8[,2] <- foreach(i=t8$letter, .combine='c') %dopar% {
- mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-</file>
-===== RgoogleMaps! =====
-The [[http://cran.r-project.org/web/packages/RgoogleMaps/index.html|RgoogleMaps]] package allows to very simply show Google maps or Google Satellite images in R, centered and zoomed on a location of your choice. You can also relatively easily overlay some spatial data from your R workspace on the map (see below). The getGeocode function tranforms a text search for a postal code or place name to latitude, longitude coordinates using Google web services.
-<file rsplus>
-library(RgoogleMaps)
-myhome=getGeoCode('Olympic stadium, Montreal');
-mymap<-GetMap(center=myhome, zoom=14)
-PlotOnStaticMap(mymap,lat=myhome['lat'],lon=myhome['lon'],cex=5,pch=10,lwd=3,col=c('red'));
-</file>
-===== Taxize =====
-The [[http://ropensci.org/|rOpenSci]] project supports the development of a number of R packages to facilitate access to a number of online data sources. Among them is the package [[http://cran.r-project.org/web/packages/taxize/index.html|Taxize]], which can be used to get taxonomic information from many different databases, including taxonomic synonyms, hierarchies, common names, and more, from a dozen different sources.
-<file rsplus>
-library(taxize)
-spp<-tax_name(query=c("american beaver"),get="species")
-fam<-tax_name(query=c("american beaver"),get="family")
-correctname <- tnrs(c("fraxinus americanus"))
-cla<-classification("acer rubrum", db = 'itis')
-</file>
-===== Spocc =====
-Another useful package from rOpenSci is [[http://cran.r-project.org/web/packages/spocc/index.html|Spocc]] which can be used to search for species occurrence data from a number of sources, including the [[http://www.gbif.org/|Global Biodiversity Information Facility]], a worldwide database containing hundreds of millions of species occurrences from collections and field data.
-<file rsplus>
-library(spocc)
-occ_data <- occ(query = 'Acer nigrum', from = 'gbif')
-mapggplot(occ_data)
-</file>
-Combine spocc and RgoogleMaps
-<file rsplus>
-occ_data <- occ(query = 'Puma concolor', from = 'gbif')
-occ_data_df=occ2df(occ_data)
-occ_data_df<-subset(occ_data_df,!is.na(latitude) & latitude!=0)
-mymap<-GetMap(center=c(mean(occ_data_df$latitude),mean(occ_data_df$longitude)), zoom=2)
-PlotOnStaticMap(mymap,lat=occ_data_df$latitude,lon=occ_data_df$longitude,cex=1,pch=16,lwd=3,col=c('red'));
-</file>
-===== geonames =====
-[[http://geonames.r-forge.r-project.org/|Geonames]] connects R to [[http://geonames.org|Geonames.org]], an online database of place names and toponyms.
-<file rsplus>
-library(geonames)
-options(geonamesUsername="glaroc")
-# Retrieve place names which contain the name "Mont Saint-Hilaire"
-res<-GNsearch(q="Mont Saint-Hilaire")
-res[,c('toponymName','fclName')]
-#Extract all citites within a bouding box defined by the coordinates of the four corners.
-dc<-GNcities(45.4, -73.55, 45.7, -73.6, lang = "en", maxRows = 10)
-dc[,c('toponymName')]
-</file>