multiknit2017.tex

%% LyX 2.0.3 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass{article}\usepackage[]{graphicx}\usepackage[]{color}
%% maxwidth is the original width if it is less than linewidth
%% otherwise use linewidth (to make sure the graphics do not exceed the margin)
\makeatletter
\def\maxwidth{ %
  \ifdim\Gin@nat@width>\linewidth
    \linewidth
  \else
    \Gin@nat@width
  \fi
}
\makeatother

\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%
\let\hlipl\hlkwb

\usepackage{framed}
\makeatletter
\newenvironment{kframe}{%
 \def\at@end@of@kframe{}%
 \ifinner\ifhmode%
  \def\at@end@of@kframe{\end{minipage}}%
  \begin{minipage}{\columnwidth}%
 \fi\fi%
 \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
 \colorbox{shadecolor}{##1}\hskip-\fboxsep
     % There is no \\@totalrightmargin, so:
     \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
 \MakeFramed {\advance\hsize-\width
   \@totalleftmargin\z@ \linewidth\hsize
   \@setminipage}}%
 {\par\unskip\endMakeFramed%
 \at@end@of@kframe}
\makeatother

\definecolor{shadecolor}{rgb}{.97, .97, .97}
\definecolor{messagecolor}{rgb}{0, 0, 0}
\definecolor{warningcolor}{rgb}{1, 0, 1}
\definecolor{errorcolor}{rgb}{1, 0, 0}
\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX

\usepackage{alltt}
\usepackage{mathpazo}
\usepackage{lmodern}
\renewcommand{\sfdefault}{lmss}
\renewcommand{\ttdefault}{lmtt}
\usepackage[T1]{fontenc}
\usepackage{geometry}
\geometry{verbose,tmargin=2.5cm,bmargin=2.5cm,lmargin=2.5cm,rmargin=2.5cm}
\usepackage{url}
\usepackage[authoryear]{natbib}

\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
% \VignetteIndexEntry{An Introduction to knitr}

\makeatother
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}
%\SweaveOpts{concordance=TRUE}


\title{Multivariate analysis in \textbf{\textsf{R}}}

\author{Oliver Purschke \\(with contributions from Helge Bruelheide, Henrik von Wehrden and Heike Zimmermann)}
\date{January 27, 2017}

\maketitle
In this practical you will perform cluster and ordination analyses, using \textbf{\textsf{R}} and the add-on packages \textbf{vegan}: \textbf{gclus}, and \textbf{clustSim}. You will use three data sets (available in \textbf{vegan}), (i) the dune meadow vegetation (\textit{dune}) and environmental data (\textit{dune.env}), (ii) a reduced version of the \textit{dune} meadow vegetation data (five of the most common species in five sites) and (iii) the \textit{varespec} and \textit{varechem} data describing vegetation and soil characteristics in a pine forest.

\section{Multivariate data}
First, lets load the required \textbf{\textsf{R}} packages and data.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(vegan)}
\hlkwd{library}\hlstd{(knitr)}
\hlkwd{library}\hlstd{(clusterSim)}
\hlkwd{library}\hlstd{(gclus)}

\hlkwd{data}\hlstd{(dune)}
\hlkwd{data}\hlstd{(dune.env)}
\end{alltt}
\end{kframe}
\end{knitrout}
\subsection{Exploring the \textit{dune} vegetation and environmental data}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{str}\hlstd{(dune)}
\end{alltt}
\begin{verbatim}
## 'data.frame':	20 obs. of  30 variables:
##  $ Achimill: num  1 3 0 0 2 2 2 0 0 4 ...
##  $ Agrostol: num  0 0 4 8 0 0 0 4 3 0 ...
##  $ Airaprae: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Alopgeni: num  0 2 7 2 0 0 0 5 3 0 ...
##  $ Anthodor: num  0 0 0 0 4 3 2 0 0 4 ...
##  $ Bellpere: num  0 3 2 2 2 0 0 0 0 2 ...
##  $ Bromhord: num  0 4 0 3 2 0 2 0 0 4 ...
##  $ Chenalbu: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Cirsarve: num  0 0 0 2 0 0 0 0 0 0 ...
##  $ Comapalu: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Eleopalu: num  0 0 0 0 0 0 0 4 0 0 ...
##  $ Elymrepe: num  4 4 4 4 4 0 0 0 6 0 ...
##  $ Empenigr: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Hyporadi: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Juncarti: num  0 0 0 0 0 0 0 4 4 0 ...
##  $ Juncbufo: num  0 0 0 0 0 0 2 0 4 0 ...
##  $ Lolipere: num  7 5 6 5 2 6 6 4 2 6 ...
##  $ Planlanc: num  0 0 0 0 5 5 5 0 0 3 ...
##  $ Poaprat : num  4 4 5 4 2 3 4 4 4 4 ...
##  $ Poatriv : num  2 7 6 5 6 4 5 4 5 4 ...
##  $ Ranuflam: num  0 0 0 0 0 0 0 2 0 0 ...
##  $ Rumeacet: num  0 0 0 0 5 6 3 0 2 0 ...
##  $ Sagiproc: num  0 0 0 5 0 0 0 2 2 0 ...
##  $ Salirepe: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Scorautu: num  0 5 2 2 3 3 3 3 2 3 ...
##  $ Trifprat: num  0 0 0 0 2 5 2 0 0 0 ...
##  $ Trifrepe: num  0 5 2 1 2 5 2 2 3 6 ...
##  $ Vicilath: num  0 0 0 0 0 0 0 0 0 1 ...
##  $ Bracruta: num  0 0 2 2 2 6 2 2 2 2 ...
##  $ Callcusp: num  0 0 0 0 0 0 0 0 0 0 ...
\end{verbatim}
\begin{alltt}
\hlkwd{str}\hlstd{(dune.env)}
\end{alltt}
\begin{verbatim}
## 'data.frame':	20 obs. of  5 variables:
##  $ A1        : num  2.8 3.5 4.3 4.2 6.3 4.3 2.8 4.2 3.7 3.3 ...
##  $ Moisture  : Ord.factor w/ 4 levels "1"<"2"<"4"<"5": 1 1 2 2 1 1 1 4 3 2 ...
##  $ Management: Factor w/ 4 levels "BF","HF","NM",..: 4 1 4 4 2 2 2 2 2 1 ...
##  $ Use       : Ord.factor w/ 3 levels "Hayfield"<"Haypastu"<..: 2 2 2 2 1 2 3 3 1 1 ...
##  $ Manure    : Ord.factor w/ 5 levels "0"<"1"<"2"<"3"<..: 5 3 5 5 3 3 4 4 2 2 ...
\end{verbatim}
\end{kframe}
\end{knitrout}
\textbf{Question:} How many species were sampled in how many sites? What character do the environmental variables have?

\subsection{Visualizing relationships among species and environmental variables}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{4}\hlstd{))}
\hlkwd{plot} \hlstd{(dune[,}\hlnum{5}\hlstd{]}\hlopt{~}\hlstd{dune[,}\hlnum{4}\hlstd{],}\hlkwc{ylab}\hlstd{=}\hlkwd{names}\hlstd{(dune[}\hlnum{5}\hlstd{]),} \hlkwc{xlab}\hlstd{=}\hlkwd{names}\hlstd{(dune[}\hlnum{4}\hlstd{]))}
\hlkwd{abline}\hlstd{(}\hlkwd{lm}\hlstd{(dune[,}\hlnum{5}\hlstd{]}\hlopt{~}\hlstd{dune[,}\hlnum{4}\hlstd{]))}
\hlkwd{plot} \hlstd{(dune[,}\hlnum{6}\hlstd{]}\hlopt{~}\hlstd{dune[,}\hlnum{4}\hlstd{],}\hlkwc{ylab}\hlstd{=}\hlkwd{names}\hlstd{(dune[}\hlnum{6}\hlstd{]),} \hlkwc{xlab}\hlstd{=}\hlkwd{names}\hlstd{(dune[}\hlnum{4}\hlstd{]))}
\hlkwd{abline}\hlstd{(}\hlkwd{lm}\hlstd{(dune[,}\hlnum{6}\hlstd{]}\hlopt{~}\hlstd{dune[,}\hlnum{4}\hlstd{]))}
\hlkwd{plot} \hlstd{(dune[,}\hlnum{5}\hlstd{]}\hlopt{~}\hlstd{dune.env[,}\hlnum{1}\hlstd{],}\hlkwc{ylab}\hlstd{=}\hlkwd{names}\hlstd{(dune[}\hlnum{5}\hlstd{]),} \hlkwc{xlab}\hlstd{=}\hlkwd{names}\hlstd{(dune.env[}\hlnum{1}\hlstd{]))}
\hlkwd{abline}\hlstd{(}\hlkwd{lm}\hlstd{(dune[,}\hlnum{5}\hlstd{]}\hlopt{~}\hlstd{dune.env[,}\hlnum{1}\hlstd{]))}
\hlkwd{plot} \hlstd{(dune[,}\hlnum{6}\hlstd{]}\hlopt{~}\hlstd{dune.env[,}\hlnum{1}\hlstd{],}\hlkwc{ylab}\hlstd{=}\hlkwd{names}\hlstd{(dune[}\hlnum{6}\hlstd{]),} \hlkwc{xlab}\hlstd{=}\hlkwd{names}\hlstd{(dune.env[}\hlnum{1}\hlstd{]))}
\hlkwd{abline}\hlstd{(}\hlkwd{lm}\hlstd{(dune[,}\hlnum{6}\hlstd{]}\hlopt{~}\hlstd{dune.env[,}\hlnum{1}\hlstd{]))}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/exploplot-1} 

\end{knitrout}
To explore many pairwise relationships between species and between species and environmental variables at once, run the following function:

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# Author: Francois Gillet, February 2007}
\hlstd{panel.cor} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{x}\hlstd{,} \hlkwc{y}\hlstd{,} \hlkwc{method} \hlstd{=} \hlstr{"pearson"}\hlstd{,} \hlkwc{digits} \hlstd{=} \hlnum{3}\hlstd{,} \hlkwc{cex.cor} \hlstd{=} \hlnum{1.2}\hlstd{) \{}
    \hlstd{usr} \hlkwb{<-} \hlkwd{par}\hlstd{(}\hlstr{"usr"}\hlstd{)}
    \hlkwd{on.exit}\hlstd{(}\hlkwd{par}\hlstd{(usr))}
    \hlkwd{par}\hlstd{(}\hlkwc{usr} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{1}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{1}\hlstd{))}
    \hlstd{r} \hlkwb{<-} \hlkwd{cor}\hlstd{(x, y,} \hlkwc{method} \hlstd{= method)}
    \hlstd{ra} \hlkwb{<-} \hlkwd{cor.test}\hlstd{(x, y,} \hlkwc{method} \hlstd{= method)}\hlopt{$}\hlstd{p.value}
    \hlstd{txt} \hlkwb{<-} \hlkwd{round}\hlstd{(r, digits)}
    \hlstd{sig} \hlkwb{<-} \hlnum{1}
    \hlstd{prefix} \hlkwb{<-} \hlstr{""}
    \hlkwa{if} \hlstd{(ra} \hlopt{<=} \hlnum{0.1}\hlstd{)}
        \hlstd{prefix} \hlkwb{<-} \hlstr{"."}
    \hlkwa{if} \hlstd{(ra} \hlopt{<=} \hlnum{0.05}\hlstd{)}
        \hlstd{prefix} \hlkwb{<-} \hlstr{"*"}
    \hlkwa{if} \hlstd{(ra} \hlopt{<=} \hlnum{0.01}\hlstd{)}
        \hlstd{prefix} \hlkwb{<-} \hlstr{"**"}
    \hlkwa{if} \hlstd{(ra} \hlopt{<=} \hlnum{0.001}\hlstd{)}
        \hlstd{prefix} \hlkwb{<-} \hlstr{"***"}
    \hlkwa{if} \hlstd{(ra} \hlopt{<=} \hlnum{0.001}\hlstd{)}
        \hlstd{sig} \hlkwb{<-} \hlnum{2}
    \hlstd{color} \hlkwb{<-} \hlnum{2}
    \hlkwa{if} \hlstd{(r} \hlopt{<} \hlnum{0}\hlstd{)}
        \hlstd{color} \hlkwb{<-} \hlnum{4}
    \hlcom{# color <- 'gray10' if(r < 0) color <- 'gray50'}
    \hlstd{txt} \hlkwb{<-} \hlkwd{paste}\hlstd{(txt, prefix,} \hlkwc{sep} \hlstd{=} \hlstr{"\textbackslash{}n"}\hlstd{)}
    \hlkwd{text}\hlstd{(}\hlnum{0.5}\hlstd{,} \hlnum{0.5}\hlstd{, txt,} \hlkwc{cex} \hlstd{= cex.cor,} \hlkwc{font} \hlstd{= sig,} \hlkwc{col} \hlstd{= color)}
\hlstd{\}}
\hlcom{## Put histograms on the diagonal}
\hlstd{panel.hist} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{x}\hlstd{,} \hlkwc{...}\hlstd{) \{}
    \hlstd{usr} \hlkwb{<-} \hlkwd{par}\hlstd{(}\hlstr{"usr"}\hlstd{)}
    \hlkwd{on.exit}\hlstd{(}\hlkwd{par}\hlstd{(usr))}
    \hlkwd{par}\hlstd{(}\hlkwc{usr} \hlstd{=} \hlkwd{c}\hlstd{(usr[}\hlnum{1}\hlopt{:}\hlnum{2}\hlstd{],} \hlnum{0}\hlstd{,} \hlnum{1.5}\hlstd{))}
    \hlstd{h} \hlkwb{<-} \hlkwd{hist}\hlstd{(x,} \hlkwc{plot} \hlstd{=} \hlnum{FALSE}\hlstd{)}
    \hlstd{breaks} \hlkwb{<-} \hlstd{h}\hlopt{$}\hlstd{breaks}
    \hlstd{nB} \hlkwb{<-} \hlkwd{length}\hlstd{(breaks)}
    \hlstd{y} \hlkwb{<-} \hlstd{h}\hlopt{$}\hlstd{counts}
    \hlstd{y} \hlkwb{<-} \hlstd{y}\hlopt{/}\hlkwd{max}\hlstd{(y)}
    \hlkwd{rect}\hlstd{(breaks[}\hlopt{-}\hlstd{nB],} \hlnum{0}\hlstd{, breaks[}\hlopt{-}\hlnum{1}\hlstd{], y,} \hlkwc{col} \hlstd{=} \hlstr{"cyan"}\hlstd{, ...)}
    \hlcom{# rect(breaks[-nB], 0, breaks[-1], y, col='gray', ...)}
\hlstd{\}}

\hlkwd{pairs}\hlstd{(}\hlkwd{cbind}\hlstd{(dune[,} \hlnum{10}\hlopt{:}\hlnum{15}\hlstd{], dune.env[,} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlnum{2}\hlstd{,} \hlnum{5}\hlstd{)]),} \hlkwc{lower.panel} \hlstd{= panel.smooth,}
    \hlkwc{upper.panel} \hlstd{= panel.cor,} \hlkwc{diag.panel} \hlstd{= panel.hist)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-3-1} 

\end{knitrout}
To quantify all possible inter-relations between the 30 species and the five environenmental descriptors, we would need to carry out 435 ($n*(n-1)/2$) separate analyses. This would be very tidious and still we wouldn't see the overall picture.

\section{Multivatiate techniques}
Multivariate methods help to  explore simultaneously the relationships between all the objects (e.g. sites) with respect to their characteristics (species composition).
\subsection{A quick appetizer}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{dune.clust} \hlkwb{<-} \hlkwd{hclust}\hlstd{(}\hlkwd{dist}\hlstd{(dune),} \hlstr{"ave"}\hlstd{)}
\hlstd{dune.ca} \hlkwb{<-} \hlkwd{decorana}\hlstd{(}\hlkwd{dist}\hlstd{(dune),} \hlkwc{ira} \hlstd{=} \hlnum{1}\hlstd{)}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{))}
\hlkwd{plot}\hlstd{(dune.clust)}
\hlkwd{plot}\hlstd{(dune.ca,} \hlkwc{type} \hlstd{=} \hlstr{"t"}\hlstd{,} \hlkwc{display} \hlstd{=} \hlstr{"sites"}\hlstd{,} \hlkwc{main} \hlstd{=} \hlstr{"Ordination"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-4-1} 

\end{knitrout}


\section{Cluster analysis (step by step)}
Here, you will perform hierarchical agglomerative clustering, which fuses objects (e.g. sites) together into single clusters that form a dendrogram.

\subsection{Example using the reduced dune data set}
Run cluster analysis on the small subset of the data (five sites x five species). Here, you use exactly the same settings (Manhattan-distance, average-linkage) that you used for generating the dendrogram in the exercise this morning.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlopt{?}\hlstd{hclust}
\hlstd{dune5} \hlkwb{<-} \hlstd{dune[}\hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{13}\hlstd{,}\hlnum{4}\hlstd{,}\hlnum{16}\hlstd{,}\hlnum{6}\hlstd{),} \hlkwd{c}\hlstd{(}\hlnum{6}\hlstd{,}\hlnum{16}\hlstd{,}\hlnum{15}\hlstd{,}\hlnum{18}\hlstd{,}\hlnum{21}\hlstd{)]}

\hlstd{dune5.man} \hlkwb{<-} \hlkwd{dist}\hlstd{(dune5,} \hlstr{"manhattan"}\hlstd{)}
\hlstd{dune5.man}
\end{alltt}
\begin{verbatim}
##     2 13  4 16
## 13  8         
## 4   1  7      
## 16  8  6  7   
## 6   8 10  7 10
\end{verbatim}
\begin{alltt}
\hlstd{dune5.man.ave.clust} \hlkwb{<-} \hlkwd{hclust}\hlstd{(dune5.man,} \hlkwc{method} \hlstd{=} \hlstr{"average"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(dune5.man.ave.clust,} \hlkwc{xlab}\hlstd{=}\hlstr{""}\hlstd{,} \hlkwc{sub}\hlstd{=}\hlstr{""}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-6-1} 

\end{knitrout}
\textbf{Q: Which two sites or groups of sites are most similar or most dissimilar? \\ Compare with the dendrogram you generated by hand. \\ Repeat the analysis with different types of linkage techniques; will those give you similar groupings of sites?}

\subsection{Using the full \textbf{dune} data set}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{dune.man} \hlkwb{<-} \hlkwd{dist}\hlstd{(dune,} \hlstr{"manhattan"}\hlstd{)}
\hlstd{dune.man.ave.clust} \hlkwb{<-} \hlkwd{hclust}\hlstd{(dune.man,} \hlkwc{method} \hlstd{=} \hlstr{"average"}\hlstd{)}
\hlkwd{plot}\hlstd{(dune.man.ave.clust,} \hlkwc{xlab}\hlstd{=}\hlstr{""}\hlstd{,} \hlkwc{sub}\hlstd{=}\hlstr{""}\hlstd{,} \hlkwc{cex}\hlstd{=}\hlnum{.8}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-7-1} 

\end{knitrout}

\subsubsection{Calculating and visualizing distance measures}
So far, you have only used the Manhattan- (City-Block) distance. Try alternative distances as well, e.g. Euclidean distance.
The \textbf{gclus} package and the coldiss - function (below) are very handy to visualize distance matrices (or square matrices in general).

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# Author: Francois Gillet, August 2009}

\hlstd{coldiss} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{D}\hlstd{,} \hlkwc{nc} \hlstd{=} \hlnum{4}\hlstd{,} \hlkwc{byrank} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{diag} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\hlstd{\{}
\hlkwd{require}\hlstd{(gclus)}
\hlkwa{if} \hlstd{(}\hlkwd{max}\hlstd{(D)}\hlopt{>}\hlnum{1}\hlstd{) D} \hlkwb{<-} \hlstd{D}\hlopt{/}\hlkwd{max}\hlstd{(D)}
\hlkwa{if} \hlstd{(byrank) \{}
\hlstd{spe.color} \hlkwb{=} \hlkwd{dmat.color}\hlstd{(}\hlnum{1}\hlopt{-}\hlstd{D,} \hlkwd{cm.colors}\hlstd{(nc))}
\hlstd{\}}
\hlkwa{else} \hlstd{\{}
\hlstd{spe.color} \hlkwb{=} \hlkwd{dmat.color}\hlstd{(}\hlnum{1}\hlopt{-}\hlstd{D,} \hlkwc{byrank}\hlstd{=}\hlnum{FALSE}\hlstd{,} \hlkwd{cm.colors}\hlstd{(nc))}
\hlstd{\}}
\hlstd{spe.o} \hlkwb{=} \hlkwd{order.single}\hlstd{(}\hlnum{1}\hlopt{-}\hlstd{D)}
\hlstd{speo.color} \hlkwb{=} \hlstd{spe.color[spe.o,spe.o]}
\hlstd{op} \hlkwb{=} \hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{),} \hlkwc{pty}\hlstd{=}\hlstr{"s"}\hlstd{)}
\hlkwa{if} \hlstd{(diag) \{}
\hlkwd{plotcolors}\hlstd{(spe.color,} \hlkwc{rlabels}\hlstd{=}\hlkwd{attributes}\hlstd{(D)}\hlopt{$}\hlstd{Labels,}
\hlkwc{main}\hlstd{=}\hlstr{"Dissimilarity Matrix"}\hlstd{,}
\hlkwc{dlabels}\hlstd{=}\hlkwd{attributes}\hlstd{(D)}\hlopt{$}\hlstd{Labels)}
\hlkwd{plotcolors}\hlstd{(speo.color,} \hlkwc{rlabels}\hlstd{=}\hlkwd{attributes}\hlstd{(D)}\hlopt{$}\hlstd{Labels[spe.o],}
\hlkwc{main}\hlstd{=}\hlstr{"Ordered Dissimilarity Matrix"}\hlstd{,}
\hlkwc{dlabels}\hlstd{=}\hlkwd{attributes}\hlstd{(D)}\hlopt{$}\hlstd{Labels[spe.o])}
\hlstd{\}}
\hlkwa{else} \hlstd{\{}
\hlkwd{plotcolors}\hlstd{(spe.color,} \hlkwc{rlabels}\hlstd{=}\hlkwd{attributes}\hlstd{(D)}\hlopt{$}\hlstd{Labels,}
\hlkwc{main}\hlstd{=}\hlstr{"Dissimilarity Matrix"}\hlstd{)}
\hlkwd{plotcolors}\hlstd{(speo.color,} \hlkwc{rlabels}\hlstd{=}\hlkwd{attributes}\hlstd{(D)}\hlopt{$}\hlstd{Labels[spe.o],}
\hlkwc{main}\hlstd{=}\hlstr{"Ordered Dissimilarity Matrix"}\hlstd{)}
\hlstd{\}}
\hlkwd{par}\hlstd{(op)}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}
Visualize Manhattan distances (magenta (reddish) - maximally similar; cyan (blueish) - maximally dissimilar):

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{coldiss}\hlstd{(dune.man,} \hlkwc{byrank} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{diag} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-9-1} 

\end{knitrout}
\textbf{Q: Compare with the Manhattan-distance matrix you calculated by hand?}
\\

Now, visualize euclidean distances:

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{dune.euclid} \hlkwb{<-} \hlkwd{dist}\hlstd{(dune,} \hlstr{"euclid"}\hlstd{)}
\hlkwd{coldiss}\hlstd{(dune.euclid,} \hlkwc{byrank} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{diag} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-10-1} 

\end{knitrout}
\textbf{Q: Are the most similar pairs of sites according to Euclidian distances similar to those obtained using Manhattan distance?}
\\
FYI: Even more distance metrics can be calculated using the "vegdist"-function in \textit{vegan}.

\subsubsection{Comparing distances - the Mantel test}
Tests correlations between two distance matrices, the extent to which one matrix resembles another one. This is somehow similar to the Pearson's correlation coeffient (r), but significance testing is done differently (using permutations).

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{mantel}\hlstd{(dune.man, dune.euclid)}
\end{alltt}
\begin{verbatim}
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = dune.man, ydis = dune.euclid) 
## 
## Mantel statistic r: 0.9514 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.105 0.144 0.174 0.214 
## Permutation: free
## Number of permutations: 999
\end{verbatim}
\end{kframe}
\end{knitrout}

\textbf{Q: Are the two distance matrices significantly similar to each other?}

\subsubsection{Finding the optimal clustering procedure}
You can use the following function to test among all possible combinations from a range of distance measures, linkage methods and number of clusters to find the optimal one.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(clusterSim)}
\hlopt{?}\hlstd{cluster.Sim}
\hlstd{clustsim} \hlkwb{<-} \hlkwd{cluster.Sim}\hlstd{(dune,} \hlnum{3}\hlstd{,} \hlnum{2}\hlstd{,} \hlnum{4}\hlstd{)}
\hlcom{# 7- mixed data (e.g. intervals, such as the cover classes in the dune dataset)}
\hlcom{# 2- minimun number of clusters, }
\hlcom{# 4 - maximum number of clusters}
\hlstd{clustsim[}\hlnum{4}\hlopt{:}\hlnum{6}\hlstd{]}
\end{alltt}
\begin{verbatim}
## $distance
##        
## "GDM2" 
## 
## $method
##           
## "average" 
## 
## $classes
## [1] "4"
\end{verbatim}
\end{kframe}
\end{knitrout}
\pagebreak
\subsection{Linking cluster analysis to environmental data}

Another clustering method (K-means clustering) can be used to generate a classification based on a prefined number of clusters (groups). Such groups can then be linked to measured environmental variables in an ANOVA-type of analysis.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{kclust} \hlkwb{<-} \hlkwd{kmeans}\hlstd{(dune,} \hlkwc{centers} \hlstd{=} \hlnum{4}\hlstd{)}

\hlkwd{par}\hlstd{(}\hlkwc{mfrow} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlnum{2}\hlstd{))}
\hlkwa{for} \hlstd{(i} \hlkwa{in} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlnum{5}\hlstd{)) \{}
    \hlkwd{boxplot}\hlstd{(dune.env[, i]} \hlopt{~} \hlstd{kclust}\hlopt{$}\hlstd{cluster,} \hlkwc{col} \hlstd{=} \hlnum{3}\hlopt{:}\hlnum{6}\hlstd{,} \hlkwc{main} \hlstd{=} \hlkwd{colnames}\hlstd{(dune.env)[i],}
        \hlkwc{xlab} \hlstd{=} \hlstr{"Cluster"}\hlstd{)}
\hlstd{\}}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-13-1} 

\end{knitrout}
\subsection{Advanced: Perform a cluster analysis on the \textit{varespec} data and link it to a continuous descriptor from the varechem data}

\pagebreak
\section{Ordination - Basic steps}
\subsection{Two species example}

Example from Legendre \& Legendre (1998), p. 392:

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{Y1} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{3}\hlstd{,}\hlnum{5}\hlstd{,}\hlnum{7}\hlstd{,}\hlnum{9}\hlstd{,}\hlnum{1}\hlstd{,}\hlnum{4}\hlstd{,}\hlnum{0}\hlstd{,}\hlnum{6}\hlstd{,}\hlnum{2}\hlstd{),}\hlkwc{nrow}\hlstd{=}\hlnum{5}\hlstd{,} \hlkwc{ncol}\hlstd{=}\hlnum{2}\hlstd{)}
\hlstd{Y1}
\end{alltt}
\begin{verbatim}
##      [,1] [,2]
## [1,]    2    1
## [2,]    3    4
## [3,]    5    0
## [4,]    7    6
## [5,]    9    2
\end{verbatim}
\begin{alltt}
\hlkwd{plot}\hlstd{(Y1,} \hlkwc{pch}\hlstd{=}\hlnum{21}\hlstd{,} \hlkwc{col}\hlstd{=}\hlstr{"black"}\hlstd{,} \hlkwc{bg}\hlstd{=}\hlstr{"black"}\hlstd{,} \hlkwc{cex}\hlstd{=}\hlnum{1.5}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-14-1} 

\end{knitrout}
\pagebreak
\subsection{Centering the attributes in a data matrix}
Shift of the origin of the coordinate system into the centre.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{Y2} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwc{nrow}\hlstd{=}\hlnum{5}\hlstd{,} \hlkwc{ncol}\hlstd{=}\hlnum{2}\hlstd{)}
\hlstd{Y2[,}\hlnum{1}\hlstd{]} \hlkwb{<-} \hlstd{Y1[,}\hlnum{1}\hlstd{]}\hlopt{-}\hlkwd{mean}\hlstd{(Y1[,}\hlnum{1}\hlstd{])}
\hlstd{Y2[,}\hlnum{2}\hlstd{]} \hlkwb{<-} \hlstd{Y1[,}\hlnum{2}\hlstd{]}\hlopt{-}\hlkwd{mean}\hlstd{(Y1[,}\hlnum{2}\hlstd{])}
\hlstd{Y2}
\end{alltt}
\begin{verbatim}
##      [,1] [,2]
## [1,] -3.2 -1.6
## [2,] -2.2  1.4
## [3,] -0.2 -2.6
## [4,]  1.8  3.4
## [5,]  3.8 -0.6
\end{verbatim}
\begin{alltt}
\hlkwd{plot}\hlstd{(Y2,} \hlkwc{pch}\hlstd{=}\hlnum{21}\hlstd{,} \hlkwc{col}\hlstd{=}\hlstr{"black"}\hlstd{,} \hlkwc{bg}\hlstd{=}\hlstr{"black"}\hlstd{,} \hlkwc{cex}\hlstd{=}\hlnum{1.5}\hlstd{,} \hlkwc{bty}\hlstd{=}\hlstr{"n"}\hlstd{,}  \hlkwc{xaxt} \hlstd{=} \hlstr{"n"}\hlstd{,} \hlkwc{yaxt} \hlstd{=} \hlstr{"n"}\hlstd{)}
\hlkwd{axis}\hlstd{(}\hlnum{1}\hlstd{,} \hlkwc{pos}\hlstd{=}\hlnum{0}\hlstd{)}
\hlkwd{axis}\hlstd{(}\hlnum{2}\hlstd{,} \hlkwc{pos}\hlstd{=}\hlnum{0}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-15-1} 

\end{knitrout}
\subsection{Rotation}
Rotate the point cloud such that the maximum variance is found on the first axis.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{dispersion.matrix} \hlkwb{<-} \hlkwd{t}\hlstd{(Y2)} \hlopt{%*%} \hlstd{Y2}

\hlcom{# divide by n-1}
\hlstd{dispersion.matrix}\hlopt{/}\hlnum{4}
\end{alltt}
\begin{verbatim}
##      [,1] [,2]
## [1,]  8.2  1.6
## [2,]  1.6  5.8
\end{verbatim}
\begin{alltt}
\hlcom{# or more simple}
\hlstd{S} \hlkwb{<-} \hlkwd{var}\hlstd{(Y2)}
\hlstd{S}
\end{alltt}
\begin{verbatim}
##      [,1] [,2]
## [1,]  8.2  1.6
## [2,]  1.6  5.8
\end{verbatim}
\begin{alltt}
\hlcom{# get the Eigenvectors}
\hlstd{U} \hlkwb{<-} \hlkwd{solve}\hlstd{(}\hlkwd{eigen}\hlstd{(S)}\hlopt{$}\hlstd{vectors)}

\hlcom{# matrix of eigenvectors U}
\hlstd{U}
\end{alltt}
\begin{verbatim}
##            [,1]       [,2]
## [1,] -0.8944272 -0.4472136
## [2,]  0.4472136 -0.8944272
\end{verbatim}
\begin{alltt}
\hlstd{F} \hlkwb{<-} \hlstd{Y2} \hlopt{%*%} \hlkwd{t}\hlstd{(U)}
\end{alltt}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(F,} \hlkwc{pch} \hlstd{=} \hlnum{21}\hlstd{,} \hlkwc{col} \hlstd{=} \hlstr{"black"}\hlstd{,} \hlkwc{bg} \hlstd{=} \hlstr{"black"}\hlstd{,} \hlkwc{cex} \hlstd{=} \hlnum{1.5}\hlstd{,} \hlkwc{bty} \hlstd{=} \hlstr{"n"}\hlstd{,} \hlkwc{xaxt} \hlstd{=} \hlstr{"n"}\hlstd{,}
    \hlkwc{yaxt} \hlstd{=} \hlstr{"n"}\hlstd{,} \hlkwc{xlim} \hlstd{=} \hlkwd{c}\hlstd{(}\hlopt{-}\hlnum{4}\hlstd{,} \hlnum{4}\hlstd{),} \hlkwc{ylim} \hlstd{=} \hlkwd{c}\hlstd{(}\hlopt{-}\hlnum{3}\hlstd{,} \hlnum{3}\hlstd{))}
\hlkwd{axis}\hlstd{(}\hlnum{1}\hlstd{,} \hlkwc{pos} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{axis}\hlstd{(}\hlnum{2}\hlstd{,} \hlkwc{pos} \hlstd{=} \hlnum{0}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-17-1} 

\end{knitrout}

\section{Unconstrained ordination (indirect gradient analysis)}
\subsection{Principal component analysis - PCA}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{data}\hlstd{(dune)}
\hlstd{model1}\hlkwb{<-}\hlkwd{rda}\hlstd{(dune,}\hlkwc{center}\hlstd{=T,}\hlkwc{scale}\hlstd{=T)}
\end{alltt}
\end{kframe}
\end{knitrout}
Remember, because it is based on Euclidean distances but species similarities among sites based on species information are not strikly Euclidean, PCA is not really suitable for species cover values. However, appropiate transformations (such as Hellinger-transformation) can circumvent this problem.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{model2} \hlkwb{<-} \hlkwd{rda}\hlstd{(}\hlkwd{decostand}\hlstd{(dune,} \hlstr{"hell"}\hlstd{),}\hlkwc{center}\hlstd{=T,}\hlkwc{scale}\hlstd{=T)}
\end{alltt}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{))}
\hlkwd{biplot}\hlstd{(model1)}
\hlkwd{biplot}\hlstd{(model2)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-20-1} 

\end{knitrout}
\textbf{Q: Does the ordination based on untransformed species data yield different configurations of sites than the ordination based on Hellinger-transformed species data?}

You can also print the summary of the eigenvalues which gives you the eigenvalues, proportion explained and cumulative proportion explained by the ordination axes (principal components):

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{summary}\hlstd{(}\hlkwd{eigenvals}\hlstd{(model2))}
\end{alltt}
\begin{verbatim}
## Importance of components:
##                         PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Eigenvalue            7.141 5.1999 3.2739 2.72312 2.42247 1.86689 1.62155
## Proportion Explained  0.238 0.1733 0.1091 0.09077 0.08075 0.06223 0.05405
## Cumulative Proportion 0.238 0.4114 0.5205 0.61127 0.69202 0.75425 0.80830
##                           PC8     PC9   PC10    PC11    PC12    PC13
## Eigenvalue            1.17624 1.03552 0.8011 0.65005 0.59471 0.53224
## Proportion Explained  0.03921 0.03452 0.0267 0.02167 0.01982 0.01774
## Cumulative Proportion 0.84751 0.88203 0.9087 0.93040 0.95022 0.96796
##                          PC14    PC15    PC16    PC17    PC18    PC19
## Eigenvalue            0.35672 0.23330 0.15034 0.11330 0.05853 0.04893
## Proportion Explained  0.01189 0.00778 0.00501 0.00378 0.00195 0.00163
## Cumulative Proportion 0.97985 0.98763 0.99264 0.99642 0.99837 1.00000
\end{verbatim}
\end{kframe}
\end{knitrout}
\textbf{Q: Interpret the summary of the eigenvectors of a PCA on the Hellinger-transformed species data.
\\
What is the percentage of total variation explained by each of first three ordination axes separately?
\\
What is the percentage of total variation explained by the first three ordination axes in conjunction?}

\subsubsection{Combining clustering and ordination results}
You can also project a dendrogram from the cluster analysis onto an ordination, to see whether sites that form a cluster in the dendrogram (i.e. are connected by blue lines) also form clusters in the ordination. And vice versa, whether sites that cluster in the ordination form a cluster in the dendrogram (i.e. are connected by blue lines).
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(model2,} \hlkwc{type} \hlstd{=} \hlstr{"p"}\hlstd{,} \hlkwc{display}\hlstd{=}\hlstr{"sites"}\hlstd{)}
\hlkwd{ordicluster}\hlstd{(model2,} \hlkwd{hclust}\hlstd{(}\hlkwd{vegdist}\hlstd{(}\hlkwd{decostand}\hlstd{(dune,} \hlstr{"hell"}\hlstd{))),} \hlkwc{prune}\hlstd{=}\hlnum{3}\hlstd{,} \hlkwc{col} \hlstd{=} \hlstr{"blue"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-22-1} 

\end{knitrout}

\textbf{Q: Are the ordination and clustering results consistent?
\\
If not, what could be the reason?}

\subsubsection{Fitting environmental vectors onto an ordination}

The interpretability of ordination plots can be enhanced by adding environmental information to an existing ordination plot.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{data}\hlstd{(dune.env)}

\hlstd{dune.env}\hlopt{$}\hlstd{n_moisture} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(dune.env}\hlopt{$}\hlstd{Moisture)}
\hlstd{dune.env}\hlopt{$}\hlstd{n_use} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(dune.env}\hlopt{$}\hlstd{Use)}
\hlstd{dune.env}\hlopt{$}\hlstd{n_manure} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(dune.env}\hlopt{$}\hlstd{Manure)}

\hlstd{fit} \hlkwb{<-} \hlkwd{envfit}\hlstd{(model2, dune.env[,}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{6}\hlstd{,}\hlnum{7}\hlstd{,}\hlnum{8}\hlstd{)],} \hlkwc{perm} \hlstd{=} \hlnum{1000}\hlstd{)}
\hlstd{fit}
\end{alltt}
\begin{verbatim}
## 
## ***VECTORS
## 
##                 PC1      PC2     r2   Pr(>r)    
## A1         -0.99604  0.08892 0.3377 0.028971 *  
## n_moisture -0.99039  0.13828 0.7860 0.000999 ***
## n_use      -0.67507 -0.73776 0.1675 0.218781    
## n_manure    0.23020 -0.97314 0.6493 0.000999 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1000
\end{verbatim}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(model2,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"sites"}\hlstd{))}
\hlkwd{plot}\hlstd{(fit)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-24-1} 

\end{knitrout}

\textbf{Q: Are the two ordination axes correlated to particular environmental gradients? \\
Identify sets of sites that have high or low values for a particular environmental variable.}

\subsubsection{PCA on environmental data}

Here, you use PCA on environmental variables to (i) reduce the multiple environmental variables into only few uncorrelated axes (dimensions) and (ii) explore correlation between environmental descriptors. Because variables are measured in different units, you need to scale the environmental variables prior to analysis.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{n_env} \hlkwb{<-} \hlstd{dune.env[,}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{6}\hlstd{,}\hlnum{7}\hlstd{,}\hlnum{8}\hlstd{)]}
\hlstd{environ} \hlkwb{<-} \hlkwd{rda}\hlstd{(n_env,}\hlkwc{center}\hlstd{=T,}\hlkwc{scale}\hlstd{=T)}
\end{alltt}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{biplot}\hlstd{(environ)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-26-1} 
\begin{kframe}\begin{alltt}
\hlkwd{summary}\hlstd{(}\hlkwd{eigenvals}\hlstd{(environ))}
\end{alltt}
\begin{verbatim}
## Importance of components:
##                         PC1    PC2    PC3     PC4
## Eigenvalue            1.684 1.3781 0.5579 0.38005
## Proportion Explained  0.421 0.3445 0.1395 0.09501
## Cumulative Proportion 0.421 0.7655 0.9050 1.00000
\end{verbatim}
\begin{alltt}
\hlkwd{scores}\hlstd{(environ,} \hlkwc{display} \hlstd{=} \hlstr{"sp"}\hlstd{)}
\end{alltt}
\begin{verbatim}
##                   PC1        PC2
## A1          1.1776484 -0.4935640
## n_moisture  1.1372860 -0.4122517
## n_use      -0.1347561 -1.3466378
## n_manure   -0.9857486 -0.8811834
## attr(,"const")
## [1] 2.952592
\end{verbatim}
\end{kframe}
\end{knitrout}

\textbf{Q: In the ordination biplot, which variables are highly correlated or uncorrelated with each other?\\
Look at the summary output. How much of the total variance in the environmental data is explained by the first two axes individually and cummulatively?\\
Which environmental gradients are represented by the first and second ordination axis, i.e. which variables have very high or very low scores on the respective axis)?}

\subsection{Correspondence analysis - CA}

Correspondence analysis has been shown be a more robust method where species show unimodal, rather than linear, responses to an underlying environmental gradient.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{CA1} \hlkwb{<-} \hlkwd{decorana}\hlstd{(dune,}\hlkwc{ira}\hlstd{=}\hlnum{1}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}
The "ira" argument defines whether standard or detrended correspondence analysis is carried out.\\
ira=1 basic reciprocal averaging = correspondence analysis= CA
ira=0, default, detrended correspondence analysis = DCA 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(CA1)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-28-1} 

\end{knitrout}
Note that site and species scores are at the same scale.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{))}
\hlkwd{plot}\hlstd{(CA1,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"sites"}\hlstd{))}
\hlkwd{plot}\hlstd{(CA1,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"species"}\hlstd{))}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-29-1} 

\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{print}\hlstd{(CA1)}
\end{alltt}
\begin{verbatim}
## 
## Call:
## decorana(veg = dune, ira = 1) 
## 
## Orthogonal correspondence analysis.
## 
##               RA1    RA2    RA3   RA4
## Eigenvalues 0.536 0.4001 0.2598 0.176
\end{verbatim}
\end{kframe}
\end{knitrout}
\textbf{Q: What are the eigenvalue for axes 1 and 2?}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{fit} \hlkwb{<-} \hlkwd{envfit}\hlstd{(CA1, n_env,} \hlkwc{perm} \hlstd{=} \hlnum{1000}\hlstd{)}
\hlstd{fit}
\end{alltt}
\begin{verbatim}
## 
## ***VECTORS
## 
##                 RA1      RA2     r2   Pr(>r)    
## A1          0.99682  0.07966 0.3104 0.035964 *  
## n_moisture  0.91712  0.39860 0.6868 0.000999 ***
## n_use       0.63528 -0.77228 0.1616 0.179820    
## n_manure   -0.12828 -0.99174 0.4848 0.000999 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1000
\end{verbatim}
\begin{alltt}
\hlkwd{scores}\hlstd{(fit,} \hlstr{"vectors"}\hlstd{)}
\end{alltt}
\begin{verbatim}
##                    RA1         RA2
## A1          0.55535366  0.04437778
## n_moisture  0.76007557  0.33034557
## n_use       0.25540159 -0.31047919
## n_manure   -0.08931823 -0.69053694
\end{verbatim}
\begin{alltt}
\hlkwd{plot}\hlstd{(CA1,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"sites"}\hlstd{))}
\hlkwd{plot}\hlstd{(fit,} \hlkwc{col} \hlstd{=} \hlstr{"blue"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-31-1} 

\end{knitrout}
\pagebreak
\subsection{Detrended correspondence analysis - DCA}

Correspondence analysis (CA) along very long environmental gradients is prone to the "arch-effect", where objects in ordination space show curved patterns. Detrended correspondence analysis (DCA) can deal with this issue.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{DCA1}\hlkwb{<-}\hlkwd{decorana}\hlstd{(dune,} \hlkwc{ira}\hlstd{=}\hlnum{0}\hlstd{,} \hlkwc{iweigh}\hlstd{=}\hlnum{1}\hlstd{)}
\hlkwd{print}\hlstd{(DCA1)}
\end{alltt}
\begin{verbatim}
## 
## Call:
## decorana(veg = dune, iweigh = 1, ira = 0) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## Downweighting of rare species from fraction 1/5.
## 
##                   DCA1   DCA2    DCA3    DCA4
## Eigenvalues     0.5033 0.2632 0.12978 0.10773
## Decorana values 0.5296 0.2243 0.04982 0.03053
## Axis lengths    3.5478 2.7827 1.48812 1.24633
\end{verbatim}
\begin{alltt}
\hlkwd{plot}\hlstd{(DCA1,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"sites"}\hlstd{),}\hlkwc{type}\hlstd{=}\hlstr{"points"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-32-1} 
\begin{kframe}\begin{alltt}
\hlstd{fit3}\hlkwb{<-} \hlkwd{envfit}\hlstd{(DCA1, n_env,} \hlkwc{perm} \hlstd{=} \hlnum{1000}\hlstd{)}
\hlkwd{plot}\hlstd{(DCA1,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"sites"}\hlstd{))}
\hlkwd{plot}\hlstd{(fit3,} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-32-2} 

\end{knitrout}
\textbf{Q: Plot CA and DCA beside each other and check whether the arch(horse-shoe)-effect disappears in DCA.}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{))}
\hlkwd{plot}\hlstd{(CA1,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"sites"}\hlstd{))}
\hlkwd{plot}\hlstd{(DCA1,,}\hlkwc{display}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"sites"}\hlstd{))}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-33-1} 

\end{knitrout}

\subsection{Principal coordinates analysis - PCoA}
While PCA assumes Euclidean distances and CA (as well as DCA) Chi-Squares distances between objects, principal coordinates analysis (PCoA) is more flexible in this respect as it allows for non-Euclidean distances.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{veg.dist} \hlkwb{<-} \hlkwd{vegdist}\hlstd{(dune,} \hlkwc{method}\hlstd{=}\hlstr{"bray"}\hlstd{)}
\hlstd{pcoa} \hlkwb{<-} \hlkwd{cmdscale}\hlstd{(veg.dist)}
\end{alltt}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(pcoa)}
\hlkwd{text}\hlstd{(pcoa[,}\hlnum{1}\hlstd{], pcoa[,}\hlnum{2}\hlstd{],} \hlkwd{rownames}\hlstd{(pcoa),} \hlkwc{cex} \hlstd{=} \hlnum{1.4}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-35-1} 

\end{knitrout}
\subsection{Nonmetric multidimensional scaling - NMDS}

Nonmetric multidimensional scaling (NMDS) is even more flexible than PCoA, as it can find non-linear relationships (unlike other methods), while being able to handle all kinds of distances.
However, you may have to run it several times to reach a good "solution". NMDS runs for a VERY long time on huge data sets.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nmds} \hlkwb{<-} \hlkwd{metaMDS}\hlstd{(dune)}
\end{alltt}
\begin{verbatim}
## Run 0 stress 0.1192678 
## Run 1 stress 0.1812938 
## Run 2 stress 0.2075713 
## Run 3 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 0.02027194  max resid 0.06496373 
## Run 4 stress 0.1192678 
## Run 5 stress 0.1192679 
## Run 6 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 2.00406e-05  max resid 5.942974e-05 
## ... Similar to previous best
## Run 7 stress 0.2003482 
## Run 8 stress 0.1183186 
## ... Procrustes: rmse 7.455498e-05  max resid 0.0001759786 
## ... Similar to previous best
## Run 9 stress 0.1183186 
## ... Procrustes: rmse 1.553252e-05  max resid 5.593033e-05 
## ... Similar to previous best
## Run 10 stress 0.1192678 
## Run 11 stress 0.1183186 
## ... Procrustes: rmse 1.626864e-05  max resid 5.213352e-05 
## ... Similar to previous best
## Run 12 stress 0.1192686 
## Run 13 stress 0.1809578 
## Run 14 stress 0.1192684 
## Run 15 stress 0.1886532 
## Run 16 stress 0.1192679 
## Run 17 stress 0.119268 
## Run 18 stress 0.1183186 
## ... Procrustes: rmse 9.306454e-05  max resid 0.0002789961 
## ... Similar to previous best
## Run 19 stress 0.1812942 
## Run 20 stress 0.1183186 
## ... Procrustes: rmse 8.355486e-05  max resid 0.0002532058 
## ... Similar to previous best
## *** Solution reached
\end{verbatim}
\begin{alltt}
\hlstd{nmds}
\end{alltt}
\begin{verbatim}
## 
## Call:
## metaMDS(comm = dune) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     dune 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1183186 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'dune'
\end{verbatim}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(nmds,} \hlkwc{type}\hlstd{=}\hlstr{"t"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-37-1} 

\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{))}
\hlkwd{plot}\hlstd{(nmds,} \hlkwc{type}\hlstd{=}\hlstr{"n"}\hlstd{)}
\hlkwd{text}\hlstd{(nmds,} \hlkwc{dis}\hlstd{=}\hlstr{"species"}\hlstd{,}\hlkwc{cex}\hlstd{=}\hlnum{0.7}\hlstd{)}
\hlcom{#points(solnmds, pch=21, col="red", bg="yellow", cex=1.2)}
\hlkwd{text}\hlstd{(nmds,} \hlstr{"sites"}\hlstd{,} \hlkwc{col}\hlstd{=}\hlstr{"blue"}\hlstd{,} \hlkwc{cex}\hlstd{=}\hlnum{0.7}\hlstd{)}
\hlkwd{plot}\hlstd{(DCA1,}\hlkwc{type}\hlstd{=}\hlstr{"n"}\hlstd{)}
\hlkwd{text}\hlstd{(DCA1,} \hlkwc{dis}\hlstd{=}\hlstr{"species"}\hlstd{,}\hlkwc{cex}\hlstd{=}\hlnum{0.7}\hlstd{)}
\hlkwd{text}\hlstd{(DCA1,} \hlkwc{dis}\hlstd{=}\hlstr{"sites"}\hlstd{,} \hlkwc{col}\hlstd{=}\hlstr{"blue"}\hlstd{,}\hlkwc{cex}\hlstd{=}\hlnum{0.7}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-38-1} 

\end{knitrout}
\subsection{Comparing ordinations using Procrustes rotation}

You can check whether two ordination are really different (i.e. have different configurations) from each other or whether differences between ordinations are just a reflection of different scalings and rotations.
Here we use the Procrustes-Rotation-Test to compare the NMDS ordination with PCA (on the Hellinger-transformed species data), and test whether the two ordinations are significantly similar.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nmds2} \hlkwb{<-} \hlkwd{metaMDS}\hlstd{(dune)}
\end{alltt}
\begin{verbatim}
## Run 0 stress 0.1192678 
## Run 1 stress 0.1192678 
## ... New best solution
## ... Procrustes: rmse 3.697152e-05  max resid 0.0001106741 
## ... Similar to previous best
## Run 2 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 0.0202829  max resid 0.06505109 
## Run 3 stress 0.1183186 
## ... Procrustes: rmse 3.742428e-05  max resid 0.0001073725 
## ... Similar to previous best
## Run 4 stress 0.1183186 
## ... Procrustes: rmse 4.915908e-06  max resid 1.398718e-05 
## ... Similar to previous best
## Run 5 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 1.981371e-05  max resid 5.996904e-05 
## ... Similar to previous best
## Run 6 stress 0.1886536 
## Run 7 stress 0.1183186 
## ... Procrustes: rmse 3.51303e-05  max resid 0.0001146371 
## ... Similar to previous best
## Run 8 stress 0.1192681 
## Run 9 stress 0.1192679 
## Run 10 stress 0.1192679 
## Run 11 stress 0.1922241 
## Run 12 stress 0.1183186 
## ... Procrustes: rmse 4.974337e-05  max resid 0.0001490783 
## ... Similar to previous best
## Run 13 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 2.01796e-05  max resid 6.444703e-05 
## ... Similar to previous best
## Run 14 stress 0.1192678 
## Run 15 stress 0.192224 
## Run 16 stress 0.1808911 
## Run 17 stress 0.1192684 
## Run 18 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 2.446088e-05  max resid 8.046263e-05 
## ... Similar to previous best
## Run 19 stress 0.1192688 
## Run 20 stress 0.1886532 
## *** Solution reached
\end{verbatim}
\begin{alltt}
\hlstd{pca2} \hlkwb{<-} \hlkwd{rda}\hlstd{(dune)}

\hlstd{protest.nmds.pca} \hlkwb{<-} \hlkwd{protest}\hlstd{(nmds2, pca2)}
\hlstd{protest.nmds.pca}
\end{alltt}
\begin{verbatim}
## 
## Call:
## protest(X = nmds2, Y = pca2) 
## 
## Procrustes Sum of Squares (m12 squared):        0.2165 
## Correlation in a symmetric Procrustes rotation: 0.8852 
## Significance:  0.001 
## 
## Permutation: free
## Number of permutations: 999
\end{verbatim}
\begin{alltt}
\hlstd{pro.nmds.pca} \hlkwb{<-} \hlkwd{procrustes}\hlstd{(nmds2, pca2)}

\hlkwd{plot}\hlstd{(pro.nmds.pca)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-39-1} 

\end{knitrout}

The following plot (Procrustes superimposition plot) shows the differences between the NMDS ordination configuration (cicles) and those that from a PCA (arrow tips).

\subsection{Surface plotting}
So far, we have only visualized linear relationships (expressed as arrows) between the objects in an ordination and environmental variables using the envfit-function, but community responses along environmental gradients are often non-linear.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{data}\hlstd{(varespec)}
\hlkwd{data}\hlstd{(varechem)}
\hlstd{nmds} \hlkwb{<-} \hlkwd{metaMDS}\hlstd{(varespec)}
\end{alltt}
\begin{verbatim}
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1843196 
## Run 1 stress 0.2290578 
## Run 2 stress 0.2114476 
## Run 3 stress 0.2048307 
## Run 4 stress 0.2218505 
## Run 5 stress 0.2251602 
## Run 6 stress 0.2414245 
## Run 7 stress 0.1825658 
## ... New best solution
## ... Procrustes: rmse 0.04169113  max resid 0.152112 
## Run 8 stress 0.2061124 
## Run 9 stress 0.195049 
## Run 10 stress 0.1845802 
## Run 11 stress 0.1967395 
## Run 12 stress 0.1969805 
## Run 13 stress 0.212657 
## Run 14 stress 0.2174195 
## Run 15 stress 0.2161251 
## Run 16 stress 0.1825664 
## ... Procrustes: rmse 0.0001286619  max resid 0.0004064493 
## ... Similar to previous best
## Run 17 stress 0.2187619 
## Run 18 stress 0.1967393 
## Run 19 stress 0.1948413 
## Run 20 stress 0.2079056 
## *** Solution reached
\end{verbatim}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{3}\hlstd{))}
\hlkwa{for} \hlstd{(i} \hlkwa{in} \hlnum{9}\hlopt{:}\hlnum{14}\hlstd{) \{}
\hlkwd{ordisurf}\hlstd{(nmds, varechem[,i],} \hlkwc{main} \hlstd{=} \hlkwd{colnames}\hlstd{(varechem)[i])}
\hlstd{\}}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-41-1} 

\end{knitrout}
\section{Constrained ordination (direct gradient analysis)}
Until now, you have performed unconstrained ordination methods, a passive form of gradient analysis, with a-posteriori interpretation of environmental variables.
Now, you will run constrained ordination, an approach that directly explores the relationships between community composition and environmental variables.


\subsection{Redundancy analysis - RDA}
RDA combines regression and PCA. Remember, when running PCA-based methods on species-data, we first have to Hellinger-transform the community data so that it represents Euclidean space.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{dune.hel} \hlkwb{<-} \hlkwd{decostand}\hlstd{(dune,} \hlstr{"hell"}\hlstd{)}
\hlstd{rda.dune} \hlkwb{<-} \hlkwd{rda}\hlstd{(dune.hel, dune.env[,}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{6}\hlopt{:}\hlnum{8}\hlstd{)])}
\hlkwd{plot}\hlstd{(rda.dune)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-42-1} 
\begin{kframe}\begin{alltt}
\hlkwd{print}\hlstd{(rda.dune)}
\end{alltt}
\begin{verbatim}
## Call: rda(X = dune.hel, Y = dune.env[, c(1, 6:8)])
## 
##               Inertia Proportion Rank
## Total          0.5559     1.0000     
## Constrained    0.2628     0.4728    4
## Unconstrained  0.2931     0.5272   15
## Inertia is variance 
## 
## Eigenvalues for constrained axes:
##    RDA1    RDA2    RDA3    RDA4 
## 0.14743 0.07987 0.02341 0.01213 
## 
## Eigenvalues for unconstrained axes:
##     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8     PC9 
## 0.06120 0.05272 0.04010 0.03481 0.02292 0.02051 0.01655 0.01324 0.00948 
##    PC10    PC11    PC12    PC13    PC14    PC15 
## 0.00701 0.00602 0.00336 0.00295 0.00142 0.00078
\end{verbatim}
\end{kframe}
\end{knitrout}

\textbf{Q: Have a look at the summary-output:\\
What is the total variance (inertia) in the data set?\\
What is the proportation of variance explained by the environmental variables (constrained proportion)?\\
What is the proportation of variance that remains unexplained (unconstrained proportion)?}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{summary}\hlstd{(}\hlkwd{eigenvals}\hlstd{(rda.dune))}
\end{alltt}
\begin{verbatim}
## Importance of components:
##                         RDA1    RDA2    RDA3    RDA4    PC1     PC2
## Eigenvalue            0.1474 0.07987 0.02341 0.01213 0.0612 0.05272
## Proportion Explained  0.2652 0.14368 0.04211 0.02181 0.1101 0.09485
## Cumulative Proportion 0.2652 0.40889 0.45099 0.47281 0.5829 0.67774
##                           PC3     PC4     PC5     PC6     PC7     PC8
## Eigenvalue            0.04010 0.03481 0.02292 0.02051 0.01655 0.01324
## Proportion Explained  0.07214 0.06261 0.04123 0.03690 0.02976 0.02381
## Cumulative Proportion 0.74988 0.81249 0.85372 0.89062 0.92039 0.94420
##                            PC9     PC10     PC11     PC12     PC13
## Eigenvalue            0.009482 0.007006 0.006018 0.003362 0.002953
## Proportion Explained  0.017060 0.012600 0.010820 0.006050 0.005310
## Cumulative Proportion 0.961250 0.973860 0.984680 0.990730 0.996040
##                           PC14      PC15
## Eigenvalue            0.001421 0.0007789
## Proportion Explained  0.002560 0.0014000
## Cumulative Proportion 0.998600 1.0000000
\end{verbatim}
\end{kframe}
\end{knitrout}
\textbf{Q: Have a look at the summary of the eigenvalues:\\
What is the proportation of variance explained by the first 2 constrained (RDA) axes?\\
What is the proportation of variance explained by the first 2 unconstrained (PC) axes for the residuals?\\
In the biplot-scores, which environmental variables are important on RDA-axis 1 and 2?\\
In the ordination-plot, which species are associated with particular management types, moisture levels etc.?
}

\pagebreak

\subsubsection{Variation partitioning}

You may also be interested in quantifying the extent to which variation in community composition is explained by different sets of variables, e.g. set 1: (soil depth (A1) and moisture (n.moisture) and set 2: land use intensity (n.use) and fertilization (manure)). Variation partitioning analysis allows partialling out the shared and unique contributions of the different sets of predictors in explaining community composition.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{soil} \hlkwb{<-} \hlstd{dune.env[,} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlnum{6}\hlstd{)]} \hlcom{# set 1}
\hlstd{management} \hlkwb{<-} \hlstd{dune.env[,} \hlkwd{c}\hlstd{(}\hlnum{7}\hlopt{:}\hlnum{8}\hlstd{)]} \hlcom{# set 2}

\hlstd{dune.part} \hlkwb{<-} \hlkwd{varpart}\hlstd{(dune.hel, soil, management)}
\hlkwd{print}\hlstd{(dune.part)}
\end{alltt}
\begin{verbatim}
## 
## Partition of variance in RDA 
## 
## Call: varpart(Y = dune.hel, X = soil, management)
## 
## Explanatory tables:
## X1:  soil
## X2:  management 
## 
## No. of explanatory tables: 2 
## Total variation (SS): 10.562 
##             Variance: 0.5559 
## No. of observations: 20 
## 
## Partition table:
##                      Df R.squared Adj.R.squared Testable
## [a+b] = X1            2   0.29705       0.21435     TRUE
## [b+c] = X2            2   0.22121       0.12958     TRUE
## [a+b+c] = X1+X2       4   0.47281       0.33222     TRUE
## Individual fractions                                    
## [a] = X1|X2           2                 0.20264     TRUE
## [b]                   0                 0.01171    FALSE
## [c] = X2|X1           2                 0.11787     TRUE
## [d] = Residuals                         0.66778    FALSE
## ---
## Use function 'rda' to test significance of fractions of interest
\end{verbatim}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(dune.part)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-45-1} 

\end{knitrout}
\textbf{Q: What is the unique contribution (individual fraction) the two variable groups (X1: soil and X2: management) in explaining community composition?\\
What percentage of the total variation remains unexplained?}

\subsection{Distance-based redundancy analysis - dbRDA}

Unlike RDA, which is based on Euclidean distances, distance-based redundancy analysis (dbRDA) can take all kinds of ecologically meaningful distances (e.g. Bray-Curtis), as it is based on principal coordinates analysis (PCoA).

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{dune.bray} \hlkwb{<-} \hlkwd{vegdist}\hlstd{(dune,} \hlstr{"bray"}\hlstd{)}
\hlstd{dune.db.rda} \hlkwb{<-} \hlkwd{capscale}\hlstd{(dune.bray} \hlopt{~} \hlstd{A1} \hlopt{+} \hlstd{n_moisture} \hlopt{+} \hlstd{n_use} \hlopt{+} \hlstd{n_manure,  dune.env)}
\end{alltt}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(dune.db.rda)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-47-1} 

\end{knitrout}
\subsection{Canonical correspondence analysis - CCA}
Canonical correspondence analysis (CCA) is the constrained version of CA, thus is based on Chi-square distances.
Unlike RDA, CCA does not allow for variation partioning analysis and doesn't have the flexibility of incorporating different types of distances (as can be done with db-RDA).
CCA-output, however, is interpretated similarily to that of RDA. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{dune.env}\hlopt{$}\hlstd{n_moisture} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(dune.env}\hlopt{$}\hlstd{Moisture)}
\hlstd{dune.env}\hlopt{$}\hlstd{n_use} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(dune.env}\hlopt{$}\hlstd{Use)}
\hlstd{dune.env}\hlopt{$}\hlstd{n_manure} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(dune.env}\hlopt{$}\hlstd{Manure)}

\hlstd{dune.cca} \hlkwb{<-} \hlstd{vegan}\hlopt{::}\hlkwd{cca}\hlstd{(dune, dune.env[,}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{6}\hlopt{:}\hlnum{8}\hlstd{)])}

\hlkwd{print}\hlstd{(dune.cca)}
\end{alltt}
\begin{verbatim}
## Call: cca(X = dune, Y = dune.env[, c(1, 6:8)])
## 
##               Inertia Proportion Rank
## Total          2.1153     1.0000     
## Constrained    0.8619     0.4075    4
## Unconstrained  1.2534     0.5925   15
## Inertia is mean squared contingency coefficient 
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4 
## 0.4411 0.2363 0.1365 0.0480 
## 
## Eigenvalues for unconstrained axes:
##     CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8     CA9 
## 0.30011 0.19798 0.17124 0.14013 0.09242 0.07666 0.06410 0.05818 0.04656 
##    CA10    CA11    CA12    CA13    CA14    CA15 
## 0.03899 0.02094 0.01708 0.01353 0.01086 0.00462
\end{verbatim}
\begin{alltt}
\hlkwd{summary}\hlstd{(}\hlkwd{eigenvals}\hlstd{(dune.cca))}
\end{alltt}
\begin{verbatim}
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CA1     CA2     CA3
## Eigenvalue            0.4411 0.2363 0.13650 0.04802 0.3001 0.19798 0.17124
## Proportion Explained  0.2085 0.1117 0.06453 0.02270 0.1419 0.09359 0.08095
## Cumulative Proportion 0.2085 0.3202 0.38475 0.40746 0.5493 0.64293 0.72388
##                           CA4     CA5     CA6     CA7     CA8     CA9
## Eigenvalue            0.14013 0.09242 0.07666 0.06410 0.05818 0.04656
## Proportion Explained  0.06625 0.04369 0.03624 0.03031 0.02750 0.02201
## Cumulative Proportion 0.79013 0.83382 0.87006 0.90036 0.92787 0.94988
##                          CA10    CA11    CA12    CA13    CA14     CA15
## Eigenvalue            0.03899 0.02094 0.01708 0.01353 0.01086 0.004618
## Proportion Explained  0.01843 0.00990 0.00807 0.00640 0.00514 0.002180
## Cumulative Proportion 0.96831 0.97821 0.98628 0.99268 0.99782 1.000000
\end{verbatim}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(dune.cca)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-49-1} 

\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{sessionInfo}\hlstd{()}
\end{alltt}
\begin{verbatim}
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.5 LTS
## 
## locale:
##  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8    
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] gclus_1.3.1       clusterSim_0.45-1 MASS_7.3-45       cluster_2.0.5    
## [5] vegan_2.4-2       lattice_0.20-34   permute_0.9-4     knitr_1.15.1     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.9     magrittr_1.5    modeest_2.1     xtable_1.8-2   
##  [5] R6_2.2.0        highr_0.6       stringr_1.1.0   tools_3.3.2    
##  [9] parallel_3.3.2  grid_3.3.2      nlme_3.1-128    mgcv_1.8-16    
## [13] e1071_1.6-7     htmltools_0.3.5 class_7.3-14    ade4_1.7-5     
## [17] digest_0.6.11   rgl_0.97.0      Matrix_1.2-8    shiny_1.0.0    
## [21] formatR_1.4     htmlwidgets_0.8 mime_0.5        evaluate_0.10  
## [25] stringi_1.1.2   R2HTML_2.3.2    jsonlite_1.2    httpuv_1.3.3
\end{verbatim}
\end{kframe}
\end{knitrout}

%\bibliographystyle{apa}
\bibliography{/home/olli/Dokumente/PhD/Literatur/Bibfile/dislit}


\end{document}