Multivariate statistics (classical methods)

Thomas Fuß

WS24/25

Multivariate statistics - in R

What you have learned so far (among other things):

A multitude of analyses how to describe the effect of one or multiple variables on one response/dependent variable


However, what can we do if we have more than one response/dependent variable?
Ecological datasets often have one or both of these datasets (if not more):

A short intro to classical ordination methods - in R

Ordination is an umbrella term for multivariate methods which arrange data along a gradient (axis, scale).

A sloppy explanation is “putting things in order”.

Ordination is a word almost exclusively used by ecologists. Outside most of these methods are known as scaling techniques.

Principal component analysis (PCA)

PCA reduces the dimensionality of a dataset while retaining most of its variability. Hence, reduces the complexity of high-dimensonal data and can identify hidden patterns.

PCA converts the original set of variables into a new set of variables called principal components (PCs), which represent the main patterns of variation.

PCA can be understood as a MLR on theoretical (latent) instead of observed dependent variables.

The regression coefficients in PCA are factor loadings.

The predicted values of the theoretical dependent are scores.

For instance, in a PCA on this matrix of n objects times p variables PC1 is computed as:

\[ PC_1=f_1X_1+f_2X_2+...+f_pX_p \] There must be as many PCs as original variables. The factor loading matrix translates p \(X\)-variables into p \(PCs\).

PCA can be understood as a rotation of the original coordinate system made up of p \(X\)-variables into a new coordinate system defined by p \(PCs\).

Principal component analysis (PCA) in R

library(vegan)
## Lade nötiges Paket: permute
## Lade nötiges Paket: lattice
## This is vegan 2.6-4
library(shape) # plots nice arrows
setwd("Z:/COURSES/datenanalyse/vu_datenanalyse/unit_3") #set working directory
mara.raw<-read.table(file="data/MaraRiver.txt",header=TRUE) # water chemistry in 54 streams, 3 types of land use
dim(mara.raw)
## [1] 54 49

head(mara.raw)
##   site sampling.time diff.solar.noon landuse order     Q canopy Temp  cond   pH
## 1 Aa20         10:04          145.92       A     2 3.572    0.2 13.4 229.1 7.45
## 2 Aa21         14:01           91.08       A     1 6.600    0.2 17.0 205.4 8.03
## 3 Aa22         16:35          245.08       A     1 1.650    0.3 14.5 205.7 7.88
## 4 Aa23         10:05          145.22       A     1 5.472    0.4 12.1 181.8 7.96
## 5 Aa6b         17:03          270.78       A     2 1.800    0.2 17.5 404.4 7.79
## 6  Af5         14:04           94.38       F     3 4.200    0.4 14.5 129.4 7.94
##     turb  Alk       TSS      PO4         NH4      NO2         NO3       Br
## 1   0.02 2160  12.55102 5.16e-07 0.000022800 8.91e-07 0.000054500 2.45e-06
## 2  35.50 2160  59.44444 8.32e-07 0.000023600 6.08e-07 0.000007810 2.06e-06
## 3   0.02 2160  10.57143 1.37e-07 0.000027800 1.13e-06 0.000047300 1.39e-06
## 4   0.02 2160  26.96629 1.00e-06 0.000020700 2.17e-06 0.000124766 1.24e-06
## 5 388.00 3740 209.41176 1.37e-07 0.000446205 1.98e-05 0.000092700 2.47e-06
## 6  28.80 2160  10.21053 1.37e-07 0.000019300 6.95e-07 0.000034900 8.64e-07
##            Cl       Fl      SO4          Na           K          Ca          Mg
## 1 0.000265196 4.47e-05 8.58e-05 0.000930796 0.000423884 0.000328085 0.000152685
## 2 0.000237582 3.91e-05 5.09e-05 0.000843062 0.000323802 0.000323170 0.000146390
## 3 0.000220969 2.87e-05 3.56e-05 0.000622227 0.000373881 0.000372199 0.000157663
## 4 0.000215891 3.03e-05 4.66e-05 0.000716442 0.000299606 0.000300339 0.000132606
## 5 0.000354103 6.12e-05 5.97e-05 0.001070465 0.000989155 0.000507910 0.000253652
## 6 0.000124757 2.87e-05 5.13e-05          NA          NA          NA          NA
##        TDN  SUVA254      slope slope_short  SR_Helms E2.to.E3 E4.to.E6      FIX
## 1 7.82e-05 2.486708 0.01484335  0.01508883 0.8385724 5.729730 7.333333 1.432674
## 2 3.20e-05 2.749288 0.01304546  0.01313918 0.8853058 4.643519 5.000000 1.352617
## 3 7.62e-05 2.957983 0.01374024  0.01376633 0.8516180 4.986425 8.000000 1.371267
## 4 1.48e-04 3.069472 0.01360719  0.01317551 0.8079401 4.892857 6.714286 1.349814
## 5 5.59e-04 1.722983 0.01444961  0.01542405 0.8943684 5.288889 7.833333 1.481364
## 6 5.49e-05 3.088685 0.01376586  0.01382645 0.8333595 4.981043 6.428571 1.336503
##       FIX2       HIX beta.alpha      X7c_c1    X7c_c2    X7c_c3    X7c_c4
## 1 1.781212 0.8178391  0.7186987 0.002585346 0.5892004 0.4082030 0.2100508
## 2 1.709486 0.7345711  0.6978052 0.075614095 0.7089239 0.4878008 0.2858983
## 3 1.726269 0.7851819  0.6829999 0.037594156 0.8459288 0.5640788 0.3114274
## 4 1.701107 0.8066407  0.6676339 0.010937179 0.6222450 0.4279429 0.2482675
## 5 1.961382 0.6234517  0.7639094 0.222123376 1.8301573 1.7719178 0.5627895
## 6 1.722147 0.8536532  0.6777914 0.035057086 1.3009351 1.0052846 0.5806384
##       X7c_c5    X7c_c6    X7c_c7 humF.to.protF humF.per.DOC protF.per.DOC
## 1 0.03701974 0.3347390 0.1175206      9.815030    0.3153769    0.03213204
## 2 0.09037746 0.3795044 0.3062165      3.943447    0.2652603    0.06726611
## 3 0.05079516 0.4395259 0.2434638      6.511800    0.3026556    0.04647802
## 4 0.04149694 0.3223648 0.1405717      8.397779    0.2706781    0.03223210
## 5 0.63429053 1.1599349 1.5369801      2.224790    0.4052359    0.18214567
## 6 0.00000000 0.7400336 0.2134428     14.595144    0.5545706    0.03799693
##    DIC.d13C     pCO2    EpCO2
## 1 -2.673634       NA       NA
## 2 -8.590664 1184.457 3.060612
## 3 -4.565283 1753.501 4.531012
## 4 -7.214370 1094.980 2.829405
## 5 -4.717236 3730.961 9.640726
## 6 -9.946630       NA       NA
mara<-mara.raw[,c(4,9:26)] #only select columns of interest

which(is.na(mara), arr.ind=TRUE) #find rows with NAs
##      row col
## [1,]   6  15
## [2,]   6  16
## [3,]   6  17
## [4,]   6  18
mara<-mara[-6,] #delete row with NAs

landuse<-mara$landuse
wc<-mara[,-1]

#apply(wc, 2, hist) #check normal distribution
wc[,-2]<-log(wc[,-2]) #log-transform concentration data (not pH)
#apply(wc, 2, hist) #check normal distribution again

plot(wc) # check potential for PCA by correlation plot




#scaling/standardizing (= subtract mean and divide by standard deviation) is necessary to account for different units
apply (wc, 2, var)
##        cond          pH        turb         Alk         TSS         PO4 
##  0.20972781  0.07256067 15.16309300  0.33465883  0.67232587  0.80413531 
##         NH4         NO2         NO3          Br          Cl          Fl 
##  0.92192840  1.08351981  1.31432386  0.46663028  0.43233487  0.14775487 
##         SO4          Na           K          Ca          Mg         TDN 
##  0.33802932  0.26713595  0.40661933  0.29809772  0.29460450  0.84457493
stand.wc <- scale(wc)
stand.wc.var <- apply (stand.wc, 2, var)
stand.wc.var
## cond   pH turb  Alk  TSS  PO4  NH4  NO2  NO3   Br   Cl   Fl  SO4   Na    K   Ca 
##    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 
##   Mg  TDN 
##    1    1

# PCA using the rda function from the vegan package
pca<-rda(X=wc,scale=TRUE)
pca<-rda(stand.wc,scale=FALSE)   # equivalent to line above

plot(pca)

biplot(pca)


scores<-scores(pca, choices=c(1:18), display="sites", scaling=1)
loadings<-scores(pca, choices=c(1:18), display="species", scaling=1) #species in our case are environmental variables

#quick sorting reveals which variables have the highest absolute correlation to the first and second axis: 
sort (abs (loadings[,1]), decreasing = TRUE)
##         Br         Na         Fl         Cl        SO4          K       cond 
## 1.80231505 1.71649988 1.66795400 1.61663710 1.58053957 1.55144854 1.48984652 
##         Mg         Ca        TDN        NO2        NH4        NO3        Alk 
## 1.48528688 1.42859134 1.29376206 1.17200054 1.15538380 1.01097565 0.94576542 
##        TSS       turb         pH        PO4 
## 0.80250922 0.78295064 0.26425447 0.08002718
sort (abs (loadings[,2]), decreasing = TRUE)
##       turb        NO3        Alk        TDN        TSS         Ca         Mg 
## 2.15423815 2.04619141 1.93776866 1.79486317 1.77397416 1.68274856 1.45997817 
##        PO4        NO2          K       cond         pH        SO4         Na 
## 1.35159505 1.33441424 1.10306586 0.79582379 0.67229923 0.55269515 0.39178510 
##         Cl         Br         Fl        NH4 
## 0.37567404 0.30952974 0.25877715 0.08609274

eigenvals(pca)
##      PC1      PC2      PC3      PC4      PC5      PC6      PC7      PC8 
## 8.061132 3.476081 1.508190 1.225212 0.715327 0.650660 0.619840 0.438682 
##      PC9     PC10     PC11     PC12     PC13     PC14     PC15     PC16 
## 0.408521 0.289131 0.166975 0.118712 0.107401 0.087748 0.063765 0.045901 
##     PC17     PC18 
## 0.010774 0.005946

# how many axes should be kept?
#there is not the one well-accepted method
#1. visual examination
#looking for a point at which the proportion of variance explained by each subsequent principal component drops off (elbow of scree plot)
screeplot(pca, npcs=length(eigenvals(pca)), type="lines") # plots eigenvalues vs. principal components (PCs) #
#2 Kaiser-Guttman criterion
# retain PCs with eigenvalues>1, PCs with eigenvalues<1 contain less information than one original variable
abline(h=1,col="red")

# eigenvalue of PC5<1
#3. variance explained method: set an "arbitrary" threshold of explained variance (e.g. 70%, 80%, 90%) and choose as many PCs till threshold is reached.
summary(pca)
## 
## Call:
## rda(X = stand.wc, scale = FALSE) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total              18          1
## Unconstrained      18          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Eigenvalue            8.0611 3.4761 1.50819 1.22521 0.71533 0.65066 0.61984
## Proportion Explained  0.4478 0.1931 0.08379 0.06807 0.03974 0.03615 0.03444
## Cumulative Proportion 0.4478 0.6410 0.72474 0.79281 0.83255 0.86870 0.90314
##                           PC8    PC9    PC10     PC11     PC12     PC13
## Eigenvalue            0.43868 0.4085 0.28913 0.166975 0.118712 0.107401
## Proportion Explained  0.02437 0.0227 0.01606 0.009276 0.006595 0.005967
## Cumulative Proportion 0.92751 0.9502 0.96627 0.975542 0.982137 0.988104
##                           PC14     PC15    PC16      PC17      PC18
## Eigenvalue            0.087748 0.063765 0.04590 0.0107740 0.0059461
## Proportion Explained  0.004875 0.003543 0.00255 0.0005986 0.0003303
## Cumulative Proportion 0.992979 0.996521 0.99907 0.9996697 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  5.531195 
## 
## 
## Species scores
## 
##           PC1      PC2      PC3      PC4      PC5      PC6
## cond -0.99702  0.34972 -0.45496 -0.04411  0.13442  0.08488
## pH    0.17684  0.29544  0.57043 -0.85920 -0.54326 -0.30724
## turb -0.52396 -0.94668 -0.34857 -0.47256 -0.03982  0.25901
## Alk  -0.63291  0.85155 -0.31425  0.13417  0.12058 -0.12839
## TSS  -0.53705 -0.77957 -0.43743 -0.64323 -0.04180  0.28207
## PO4  -0.05355 -0.59396  0.73584 -0.33539  0.79636 -0.11188
## NH4  -0.77319  0.03783 -0.68069 -0.06619  0.19035 -0.53093
## NO2  -0.78431 -0.58641 -0.15264 -0.08530 -0.10877 -0.62400
## NO3  -0.67655 -0.89920  0.22318  0.46162 -0.18975 -0.01087
## Br   -1.20613 -0.13602  0.08215  0.06930  0.01147  0.28837
## Cl   -1.08187 -0.16509  0.49846  0.21967  0.03718 -0.12252
## Fl   -1.11621 -0.11372  0.18254 -0.05032  0.14197 -0.03546
## SO4  -1.05771  0.24288  0.40207  0.23188 -0.21474  0.01286
## Na   -1.14870  0.17217  0.13648  0.05892 -0.21020  0.18551
## K    -1.03824  0.48474  0.25134 -0.15451  0.03470  0.09851
## Ca   -0.95603  0.73948 -0.06107 -0.22861  0.11002  0.06037
## Mg   -0.99397  0.64159  0.05003 -0.24907  0.04694  0.06880
## TDN  -0.86580 -0.78875 -0.05764  0.31464 -0.19754 -0.07204
## 
## 
## Site scores (weighted sums of species scores)
## 
##          PC1      PC2       PC3      PC4      PC5      PC6
## 1  -1.283073  1.46362  0.329969  1.17010  0.42054  0.64811
## 2  -1.034856  1.29806  0.047297 -1.40897  0.23982  1.05835
## 3  -0.882710  1.69037 -0.205471  0.53435 -0.97579 -0.24031
## 4  -0.994400  0.94315  0.505607  0.14215 -0.25215 -0.60337
## 5  -2.431354  0.65908 -1.898250 -0.84073 -0.96461 -1.22401
## 7   0.625497  0.59533  0.335832  0.26999  0.46695 -1.53700
## 8   1.479591  0.71474 -2.670917 -1.30598 -0.33714  0.79142
## 9  -0.445864  0.46757 -0.662112  2.69914 -0.97712  0.59770
## 10 -0.761641  0.42997  0.078779 -1.25501  1.05851  0.23628
## 11  0.236561  0.72111 -0.598763  0.65145  1.19091  0.35206
## 12  1.141965  0.47600 -0.093490  0.25066 -0.74183 -0.50609
## 13 -0.266793 -0.81652 -0.519182 -0.08502  0.98419  0.15946
## 14  1.856240 -1.27798 -2.134725  0.31427  0.53622 -1.01871
## 15 -0.551353 -1.32214 -0.026774 -0.19529 -0.31291  0.30059
## 16 -1.000529 -0.09745 -0.165024  0.82015  0.76709 -1.46104
## 17 -1.137819 -0.55812  0.102472  0.12211  0.28418 -0.34966
## 18 -0.635730 -0.10995 -0.083433  0.08852  0.08615  1.19206
## 19 -0.350048 -0.61812  0.238228 -0.74534 -0.27547  0.09870
## 20  0.043308 -0.44619 -0.616293  0.21443 -1.00207  1.03184
## 21 -0.506526 -0.74355  0.520934 -0.32837 -0.52971 -0.12075
## 22 -0.191623 -0.73919  0.679868 -0.07152 -0.12217  2.02785
## 23 -0.542921 -0.42120  0.205271 -0.34908  0.08852 -0.25500
## 24 -0.891549 -1.03600  0.087247 -0.06393  0.27935 -0.67049
## 25 -0.046249 -0.83270  0.886679 -1.66658 -1.29329 -1.16661
## 26  0.287167 -1.00938  0.075905 -0.38458  0.26442  0.03373
## 27  0.062323 -0.58078 -0.420579  0.55599 -0.74345  0.20414
## 28  0.114344 -1.02675 -0.128087 -0.11522  0.48227 -0.18664
## 29  0.131909 -0.90355  0.092727  0.14660  0.57868  0.27300
## 30 -0.242136 -0.88933 -0.900465  0.49281 -0.20915  0.35172
## 31  0.598040  0.66619  0.632611  0.63693 -0.11075  1.10934
## 32  0.291300  0.63471  0.038214  0.41094  0.04070 -0.21870
## 33  0.143544 -0.10039  0.902594  1.70768 -0.31291 -1.76226
## 34 -0.123959  0.49543  0.498420 -1.42250  0.43705 -0.16622
## 35 -0.007017 -0.44478  0.598281 -0.30532  0.06941 -0.01775
## 36  0.325450  0.81900 -0.007016 -0.05172 -0.22195 -0.04243
## 37  0.792494  0.37487 -0.251334  0.32512 -0.90002  0.48860
## 38  0.578051  0.34156  0.563571 -0.15707  0.78511 -0.50635
## 39  0.678964  0.74658  0.548865  0.21253  1.15356  0.84082
## 40  0.552941  0.42893  0.518675 -0.06526  1.10748 -0.31764
## 41  1.096463  0.22508  1.878217 -0.98915 -1.08047 -0.53829
## 42  1.122814  0.65436 -0.062654  0.16617 -2.25180 -0.79297
## 43  0.388900  0.62223  0.979109 -0.59038  0.67362  0.28388
## 44  0.936463  0.08358  0.377874  0.30759  0.86156 -0.93175
## 45  0.630656  1.52392 -1.054410 -1.08937  0.97352 -0.21177
## 46 -0.460918 -0.86694 -0.540394 -0.12777  0.13635 -0.42237
## 47 -0.135065 -0.19959 -0.178417 -0.17257  0.24690 -0.55419
## 48  0.363278  0.39923  0.192034  0.62466  0.65023 -0.15452
## 49 -0.146583 -0.07110 -0.609688 -0.31289 -0.43489  0.59717
## 50 -0.184906 -0.65088 -0.098466  0.17505 -0.44824  0.50493
## 51 -0.329674 -0.27009  0.717715  0.25504 -0.12641  0.24305
## 52  0.668426 -0.10346  0.834806 -0.37667 -1.70236  1.04501
## 53  0.230527 -0.56106  0.072573  0.51022  0.87297  0.49728
## 54  0.208081 -0.77747  0.385568  0.67163  0.59043  1.00980
# -> the first 3-4 PCs seem useful, and just PC1 and PC2 alone are already explaining a lot of overall variance


#Different scaling produces different plots
#Either species/environmental variables (2) or site (1) scores are scaled by eigenvalues, and the other set of scores is left unscaled, or with 3 both are scaled symmetrically by square root of eigenvalues.
# for a DISTANCE BIPLOT (focus is on sites, "scaling 1")
# in this plot:
# 1) distances among sites are approximating true Euclidean distances in multivariate space
# 2) angles between arrows do not reflect correlations among variables
# 3) projecting site on descriptor at right angle gives its appr. descriptor value
# each principal component has variance given by eigenvalue, loadings remain unscaled
str(landuse)
##  chr [1:53] "A" "A" "A" "A" "A" "M" "F" "M" "M" "M" "F" "A" "A" "A" "A" "A" ...
landuse<-as.factor(landuse)

plot(scores[,1:2],asp=1,pch=21, cex=2, bg=landuse) # note asp=1
arrows<-loadings[,1:2]*0.7 # with extension/reduction factor
Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="darkgreen")
text(x=arrows[,1]*1.2,y=arrows[,2]*1.2,labels=names(wc),cex=0.8)
#text(scores[,1:2], rownames(wc), pos=2)
legend("bottomright", legend=c("Agriculture", "Mixed", "Forest"), pch=21, pt.cex=1, pt.bg=c("black","green","darkred"), cex=0.8)


#biplot(pca,scaling=1)

# for a CORRELATION BIPLOT (focus is on variables, "scaling 2")
# in this plot
# 1) distances among sites are not approximating true Euclidean distances in multivariate space
# 2) angles between arrows reflect correlations among variables (NOT proximity of arrow heads)
# 3) projecting site on descriptor at right angle gives its appr. descriptor value
# each principal component is weighted by 1/sqrt(eigenvalue), so it has variance 1
scores<-scores(pca, choices=c(1:18), display="sites", scaling=2)
loadings<-scores(pca, choices=c(1:18), display="species", scaling=2) #species in our case are environmental variables

plot(scores[,1:2],asp=1,pch=21, cex=2, bg=landuse) # note asp=1
arrows<-loadings[,1:2] 
Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="darkgreen")
text(x=arrows[,1]*1.2,y=arrows[,2]*1.2,labels=names(wc),cex=0.8)
#text(scores[,1:2], rownames(wc), pos=2)
legend("bottomright", legend=c("Agriculture", "Mixed", "Forest"), pch=21, pt.cex=1, pt.bg=c("black","green","darkred"), cex=0.8)


#biplot(pca,scaling=2)#distances among objects are not approximations of Euclidean distances; angles between descriptor (species) vectors reflect their correlations.
# Descriptors at 180 degrees of each other are negatively correlated;
# Descriptors at 90 degrees of each other have zero correlation;
# Descriptors at 0 degrees of each other are positively correlated.



#Scaling 3 is a compromise
scores<-scores(pca, choices=c(1:18), display="sites", scaling=3)
loadings<-scores(pca, choices=c(1:18), display="species", scaling=3) #species in our case are environmental variables

plot(scores[,1:2],asp=1,pch=21, cex=2, bg=landuse) # note asp=1
arrows<-loadings[,1:2] 
Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="darkgreen")
text(x=arrows[,1]*1.2,y=arrows[,2]*1.2,labels=names(wc),cex=0.8)
#text(scores[,1:2], rownames(wc), pos=2)
legend("bottomright", legend=c("Agriculture", "Mixed", "Forest"), pch=21, pt.cex=1, pt.bg=c("black","green","darkred"), cex=0.8)


#biplot(pca,scaling=3)


# Alternatively:
# #Run PCA with R basic
# pca<-prcomp(scale(wc),center=F,scale.=F) #scaling (= subtract mean and divide by standard deviation) is necessary to account for different units
# pca<-prcomp(wc,center=T,scale.=T)   # equivalent to line above
# 
# summary(pca)
# str(pca)
# 
# pca$sdev # stdevs of PCs (squares are eigenvalues)
# eigenvals(pca)==pca$sdev^2 #eigenvalues correspond to the variance (=sd^2) of their respective principal component (PC)
# sum(eigenvals(pca)) #total variance of the dataset
# 
# head(scores<-pca$x) # site scores on all PCs
# #head(scale(wc) %*% pca$rotation) # to manually compute scores from variables and loadings
# 
# head(loadings<-pca$rotation) # the variable loadings
# 
# 
# # for a DISTANCE BIPLOT (focus is on sites, "scaling 1")
# # in this plot:
# # 1) distances among sites are approximating true Euclidean distances in multivariate space
# # 2) angles between arrows do not reflect correlations among variables
# # 3) projecting site on descriptor at right angle gives its appr. descriptor value
# # each principal component has variance given by eigenvalue, loadings remain unscaled
# str(landuse)
# landuse<-as.factor(landuse)
# 
# layout(matrix(1:2, nrow=1))
# plot(scores[,1:2],asp=1,pch=21, cex=2, bg=landuse) # note asp=1
# arrows<-loadings*15 # with extension factor
# Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="darkgreen")
# text(x=arrows[,1]*1.2,y=arrows[,2]*1.2,labels=names(wc),cex=0.8)
# #text(scores[,1:2], rownames(wc), pos=2)
# legend("bottomright", legend=c("Agriculture", "Mixed", "Forest"), pch=21, pt.cex=1, pt.bg=c("black","green","darkred"), cex=0.8)
# 
# biplot(pca,scale=0)
# 
# # for a CORRELATION BIPLOT (focus is on variables, "scaling 2")
# # in this plot
# # 1) distances among sites are not approximating true Euclidean distances in multivariate space
# # 2) angles between arrows reflect correlations among variables (NOT proximity of arrow heads)
# # 3) projecting site on descriptor at right angle gives its appr. descriptor value
# # each principal component is weighted by 1/sqrt(eigenvalue), so it has variance 1
# var(scores[,1]/pca$sdev[1]) # just demo
# plot(scores[,1]/pca$sdev[1],scores[,2]/pca$sdev[2],pch=21,bg=landuse,asp=1)
# # loadings are weighted by sqrt(eigenvalues)
# arrows<-loadings*matrix(pca$sdev,nrow=nrow(loadings),ncol=ncol(loadings),byrow=TRUE)
# arrows<-arrows*2 # choose extension factor
# Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="purple")
# # as alternative just compute correlation of scores with original data ("structure coefficients")
# (structure<-cor(wc,scores))
# structure<-2*structure
# Arrows(x0=0,y0=0,x1=structure[,1],y1=structure[,2],col="red")
# text(x=arrows[,1]*1.3,y=arrows[,2]*1.2,labels=names(wc),cex=0.7)
# 
# biplot(pca,scale=1)


##############################
# some follow-up suggestions #
# test PCA-axes for effect of landuse or stream size (as log(Q)) using ANOVA or ANCOVA
# correlate PCA-axes with other potential "controlling" variables (e.g. canopy cover) to give "meta-dimensions" more meaning
# useful function envfit() to relate additional variables to the ordination space

# The function envfit calculates multiple regression of environmental variable with ordination axes (environmental variable is used as dependent and selected ordination axes as explanatory variables). Significance is tested by permutation test. Vectors (for continual variables) and centroids (for categorical variables) can be projected onto ordination diagram using plot function.

Direct (constrained) ordination

The classical ordination methods always target dimension reduction and may involve indirect (unconstrained) reconstruction of gradients from a single data matrix (e.g. PCA) or direct (constrained) regression on a 2nd matrix of explanatory variables (e.g. RDA).

Constrained ordinations aim to explain the variation in a set of response variables using a set of explanatory variables.


Redundancy analysis (RDA)

Direct ordinations identify the relationship between a multivariate response matrix and a multivariate explanatory matrix by combining regressino and ordination concepts.

Two involved matrices: one dependent, one independent.

Steps in RDA

  1. Multiple linear regression relate each dependent/response variable (eg. species abundances, DOM descriptors) to the independent/explanatory matrix (e.g. environmental variables, land use) and predict the response.

  2. The matrix of predicted response variables (same size as original: number of sites/samples/objects * number of variables) is used in PCA to extract constrained ordination axes. Hence, RDA is a canonical version of PCA where the PCs are constrained to be linear combinations of the explanatory variables.

  3. The matrix of the residuals is also used in PCA to extract unconstrained axes.

Differences to PCA:

  1. The RDA can only “ordinate” variation of the responses that is relatable to predictors. Hence, gradients in the response matrix, which cannot be related to the explanatory matrix (as linear combinations), cannot be identified. Assume the most important environmental variable was not measured in the sampling campaign -> compare PCA to RDA.

  2. Often we are interested in the amount and significance of variation in the response matrix (e.g. species composition) that can be explained by the explanatory matrix (e.g. environmental variables). Significance tests are available (permutation-based).

For simplicity, the green (=independent) matrix only contains one variable, normally you have more than one

Redundancy analysis (RDA) in R

zwc<-scale(wc) # must at least be centered even if dimensionally homogeneous!
xmat<-data.frame(logQ=log(mara.raw$Q),logTDN=log(mara.raw$TDN),canopy=mara.raw$canopy)[-6,]
xmat<-data.frame(scale(xmat))

rda<-rda(zwc~logQ+logTDN+canopy,data=xmat) # take care: confusing X and Y argument names

# actual RDA output check
summary(rda)
## 
## Call:
## rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          18.000       1.00
## Constrained     6.481       0.36
## Unconstrained  11.519       0.64
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                         RDA1    RDA2    RDA3    PC1     PC2     PC3     PC4
## Eigenvalue            5.5244 0.76796 0.18831 4.9837 1.63026 1.19739 0.77385
## Proportion Explained  0.3069 0.04266 0.01046 0.2769 0.09057 0.06652 0.04299
## Cumulative Proportion 0.3069 0.34957 0.36004 0.6369 0.72748 0.79400 0.83699
##                           PC5     PC6     PC7     PC8     PC9     PC10     PC11
## Eigenvalue            0.64230 0.60556 0.46542 0.38373 0.25602 0.153021 0.117674
## Proportion Explained  0.03568 0.03364 0.02586 0.02132 0.01422 0.008501 0.006537
## Cumulative Proportion 0.87268 0.90632 0.93217 0.95349 0.96772 0.976216 0.982754
##                           PC12     PC13     PC14     PC15      PC16      PC17
## Eigenvalue            0.101408 0.079286 0.061438 0.046321 0.0156525 0.0063286
## Proportion Explained  0.005634 0.004405 0.003413 0.002573 0.0008696 0.0003516
## Cumulative Proportion 0.988387 0.992792 0.996205 0.998779 0.9996484 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1   RDA2    RDA3
## Eigenvalue            5.5244 0.7680 0.18831
## Proportion Explained  0.8524 0.1185 0.02906
## Cumulative Proportion 0.8524 0.9709 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  5.531195 
## 
## 
## Species scores
## 
##         RDA1      RDA2      RDA3        PC1        PC2        PC3
## cond -0.6632  0.526208 -0.053197 -7.126e-01 -2.642e-01 -2.366e-01
## pH    0.4359 -0.080604 -0.110311 -2.354e-01  6.035e-01 -2.708e-01
## turb -0.8597  0.004151  0.171732  3.490e-01  4.901e-01 -5.836e-01
## Alk  -0.1085  0.409957 -0.026803 -8.856e-01 -5.159e-01 -4.483e-02
## TSS  -0.7648  0.195626  0.067225  2.585e-01  4.908e-01 -7.306e-01
## PO4  -0.1714 -0.057514 -0.327410  2.628e-01  9.564e-01  3.017e-01
## NH4  -0.6307  0.217307 -0.051937 -4.163e-01 -3.196e-01 -5.981e-01
## NO2  -0.8571 -0.178053  0.281840 -1.528e-01  2.746e-01 -3.129e-01
## NO3  -1.1170 -0.494302 -0.147602  1.649e-01  1.328e-01  3.170e-01
## Br   -0.9585  0.116087 -0.072049 -6.810e-01  2.180e-01  6.036e-02
## Cl   -0.8351 -0.016235  0.006982 -6.587e-01  4.095e-01  4.346e-01
## Fl   -0.7748  0.116617 -0.017050 -7.011e-01  4.096e-01  3.659e-02
## SO4  -0.5605  0.008991  0.140803 -9.589e-01  2.172e-01  3.459e-01
## Na   -0.7506  0.071258  0.039663 -8.902e-01  1.087e-01  6.183e-02
## K    -0.4348  0.085161 -0.092075 -1.085e+00  1.745e-01 -4.157e-02
## Ca   -0.3188  0.430074 -0.088844 -1.078e+00 -7.507e-02 -1.775e-01
## Mg   -0.3587  0.282375 -0.044108 -1.095e+00  4.243e-02 -1.784e-01
## TDN  -1.2252 -0.429552 -0.118224  9.845e-17 -1.289e-16  1.410e-16
## 
## 
## Site scores (weighted sums of species scores)
## 
##         RDA1      RDA2     RDA3      PC1      PC2       PC3
## 1  -0.899757  3.084889 -0.55294 -1.73973 -0.49579  1.405403
## 2  -0.601643  4.018791  0.23349 -1.62902  1.01934 -0.493815
## 3  -0.425217  2.834368  0.10401 -1.40811 -1.10996  0.502236
## 4  -0.756630  1.881205 -1.51654 -1.37687 -0.19517  0.362439
## 5  -2.455504  4.254210  2.42612 -2.05980 -0.84349 -2.485533
## 7   0.861601 -0.365693 -1.72079 -0.14991 -0.56146  0.155296
## 8   1.764345  1.159717  2.64686  0.62914 -1.54289 -2.968104
## 9  -0.482278  0.089585  0.15601 -0.35169 -2.59873  1.153105
## 10 -0.595642  2.398130 -0.19763 -0.94864  1.32755 -0.723328
## 11  0.451013  1.306409 -1.18506 -0.42811 -0.88023 -0.088124
## 12  1.394450 -0.925289  1.78413 -0.03078 -0.38908 -0.245790
## 13 -0.567171 -0.078456  0.53359  0.88444  0.28402  0.167429
## 14  1.469933 -2.003658  3.48847  3.00526 -1.06668 -0.266833
## 15 -1.033131 -1.201150  0.27119  0.59853  0.67023 -0.195910
## 16 -1.140389  0.708121  0.03320 -0.37468 -0.34163  0.496573
## 17 -1.426847  0.273804 -0.49117 -0.56728  0.24291 -0.207985
## 18 -0.726797  0.919409 -0.63159 -0.09927  0.15910  0.396983
## 19 -0.562127 -0.441430 -0.39299  0.45991  0.63986 -0.006655
## 20 -0.146589 -0.580216  1.46528  0.62302 -0.45546 -0.030056
## 21 -0.792760 -0.936524 -1.52046  0.06731  0.43021 -0.279145
## 22 -0.450183 -0.920766 -1.69925  0.35532  0.44301  0.394593
## 23 -0.722434 -0.099725 -1.05913  0.11364  0.32640  0.053559
## 24 -1.294677 -0.621767  0.54340 -0.14052  0.78301 -0.387924
## 25 -0.256636 -1.469854 -1.00794  0.30718  1.42878 -0.901286
## 26 -0.006011 -1.524849  0.44526  0.71389  0.74750 -0.370880
## 27 -0.167762 -1.210237  1.92863  0.73162 -0.47637  0.287856
## 28 -0.214533 -1.385767  0.40575  0.68065  0.48451 -0.293713
## 29 -0.153241 -1.359169  0.41669  0.54181  0.45796  0.113385
## 30 -0.608704 -0.767504  2.40998  0.58176 -0.34287 -0.271359
## 31  0.877051  0.127594 -0.99204 -0.11448 -0.08793  1.149609
## 32  0.524430  0.342174  0.11586 -0.19013 -0.41330  0.388517
## 33  0.100447 -1.912077 -0.27986 -0.07708 -0.63553  1.401936
## 34  0.113385  1.339155 -0.76471 -0.15951  1.13570 -0.158967
## 35 -0.144093 -0.913297 -2.68176  0.14873  0.34036 -0.116528
## 36  0.635959  0.706963 -0.20569 -0.45058 -0.27741 -0.210309
## 37  0.964039 -0.502987  1.56066  0.58492 -0.58515  0.246218
## 38  0.780714 -0.223068 -1.48202 -0.02480  0.36112  0.196388
## 39  1.002097  0.320595 -2.44794 -0.11533 -0.06812  0.544981
## 40  0.777442 -0.009334 -1.95549 -0.09691  0.23264  0.233320
## 41  1.374032 -1.702035  0.01784 -0.03704  1.50770  0.434674
## 42  1.421599 -1.106980  2.45558 -0.07522 -0.69935 -0.519837
## 43  0.700872  0.512436 -1.78976 -0.23677  0.87534  0.442824
## 44  1.074898 -1.074960  0.76215  0.24269  0.26216  0.475054
## 45  1.160824  2.438268 -2.38973  0.14196 -1.14818 -0.702500
## 46 -0.808194 -0.336965  0.64062 -0.06377  0.17315 -1.013755
## 47 -0.197930  0.127281  0.76092 -0.11609  0.31741 -0.241202
## 48  0.528880 -0.060544  0.32122  0.08513 -0.24501  0.789760
## 49 -0.190294  0.489138  1.74631 -0.11103  0.04844 -0.769059
## 50 -0.443365 -0.924883 -0.08718  0.49389 -0.16078  0.116308
## 51 -0.445001 -0.638721 -2.14463  0.13713  0.06390  0.808019
## 52  0.728592 -1.401673  2.01083  0.17276  0.63744  0.369228
## 53  0.045616 -1.090325 -1.06195  0.41535 -0.04511  0.260961
## 54 -0.036682 -1.542340  0.57423  0.45711  0.26592  0.601943
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##         RDA1     RDA2     RDA3      PC1      PC2       PC3
## 1  -0.192647  0.93038  0.92799 -1.73973 -0.49579  1.405403
## 2   0.427478  1.34335  1.22481 -1.62902  1.01934 -0.493815
## 3  -0.146209  0.99448  0.45228 -1.40811 -1.10996  0.502236
## 4  -0.310322 -0.13504  0.14698 -1.37687 -0.19517  0.362439
## 5  -1.484641 -0.21700  0.38344 -2.05980 -0.84349 -2.485533
## 7   0.674710 -0.78138 -0.51445 -0.14991 -0.56146  0.155296
## 8   1.388710 -0.10748 -0.22376  0.62914 -1.54289 -2.968104
## 9  -0.858711 -0.86050 -1.00156 -0.35169 -2.59873  1.153105
## 10  0.407941  1.11874 -0.68765 -0.94864  1.32755 -0.723328
## 11  0.703912 -0.23995 -0.63485 -0.42811 -0.88023 -0.088124
## 12  1.393070 -1.09377  1.50118 -0.03078 -0.38908 -0.245790
## 13 -0.864232  1.33162 -0.10131  0.88444  0.28402  0.167429
## 14  0.070930  1.59769  0.43814  3.00526 -1.06668 -0.266833
## 15 -1.113644 -0.02139  0.08487  0.59853  0.67023 -0.195910
## 16 -0.972647  0.47906 -0.44170 -0.37468 -0.34163  0.496573
## 17 -1.040565 -0.76624 -0.16526 -0.56728  0.24291 -0.207985
## 18 -0.629396  0.77943  0.19352 -0.09927  0.15910  0.396983
## 19 -0.880788  0.52142  0.58591  0.45991  0.63986 -0.006655
## 20 -0.557724  0.29157  0.85510  0.62302 -0.45546 -0.030056
## 21 -0.816524 -0.74571 -1.49309  0.06731  0.43021 -0.279145
## 22 -0.679726 -0.61138 -0.94934  0.35532  0.44301  0.394593
## 23 -0.876393  0.22614  0.17538  0.11364  0.32640  0.053559
## 24 -0.835658 -0.46318 -0.11007 -0.14052  0.78301 -0.387924
## 25 -0.306131 -0.56840 -0.71753  0.30718  1.42878 -0.901286
## 26 -0.050644 -0.33280 -0.12269  0.71389  0.74750 -0.370880
## 27 -0.663953  0.25632  0.79512  0.73162 -0.47637  0.287856
## 28 -0.304155 -0.31610  0.19760  0.68065  0.48451 -0.293713
## 29 -0.229915 -0.50523  0.78395  0.54181  0.45796  0.113385
## 30 -0.750822  0.01951  0.79961  0.58176 -0.34287 -0.271359
## 31  0.919784  0.24149 -0.12906 -0.11448 -0.08793  1.149609
## 32  0.500061  0.28682  0.07581 -0.19013 -0.41330  0.388517
## 33 -0.169995 -1.32278  0.05458 -0.07708 -0.63553  1.401936
## 34  0.327650  1.71656 -0.89375 -0.15951  1.13570 -0.158967
## 35 -0.278900 -0.83339 -1.12100  0.14873  0.34036 -0.116528
## 36  0.827677  0.04073 -0.62665 -0.45058 -0.27741 -0.210309
## 37  0.573226  0.88396 -1.00915  0.58492 -0.58515  0.246218
## 38  0.850153 -0.05214 -0.58813 -0.02480  0.36112  0.196388
## 39  1.044336  0.18841 -1.50261 -0.11533 -0.06812  0.544981
## 40  0.840242 -0.12799 -0.57504 -0.09691  0.23264  0.233320
## 41  1.451280 -0.87165  0.50304 -0.03704  1.50770  0.434674
## 42  1.324854 -1.06452 -0.01889 -0.07522 -0.69935 -0.519837
## 43  0.954237  0.72568 -1.21010 -0.23677  0.87534  0.442824
## 44  1.025094 -0.47834  1.59291  0.24269  0.26216  0.475054
## 45  0.465138  1.60579  0.20294  0.14196 -1.14818 -0.702500
## 46 -0.461071 -1.14745 -0.17978 -0.06377  0.17315 -1.013755
## 47  0.002323 -0.19579  1.34466 -0.11609  0.31741 -0.241202
## 48  0.426895  0.53629 -0.03598  0.08513 -0.24501  0.789760
## 49  0.119534  0.16802 -0.14099 -0.11103  0.04844 -0.769059
## 50 -0.832201 -0.18606  0.31023  0.49389 -0.16078  0.116308
## 51 -0.803215 -0.10485 -0.18123  0.13713  0.06390  0.808019
## 52  0.580101 -0.79935  1.87566  0.17276  0.63744  0.369228
## 53 -0.091596 -0.62895 -0.07307  0.41535 -0.04511  0.260961
## 54 -0.096912 -0.70466 -0.05704  0.45711  0.26592  0.601943
## 
## 
## Biplot scores for constraining variables
## 
##           RDA1    RDA2     RDA3 PC1 PC2 PC3
## logQ    0.5174 -0.6315  0.57757   0   0   0
## logTDN -0.9398 -0.3295 -0.09068   0   0   0
## canopy  0.5328 -0.2110 -0.81953   0   0   0


#Break it down:
# Constrained Proportion: variance of Y explained by X (36%)
sum(eigenvals(rda)[1:3])/sum(eigenvals(rda))
## [1] 0.3600362
# Unconstrained Proportion: unexplained variance in Y (64%)

# How would you report these results? You could say: “The included predictors explain 36% of the variation in water chemistry across sites.”
# 
# Summary also contains:
# Eigenvalues, and their contribution to the variance
# Scores for species, sites, and the explanatory variables, which are the coordinates of each of these objects in the RDA space. The default scaling is of type 2.
# 
# 

RsquareAdj(rda) # Adjusted R2 measures the strength of the relationship between Y and X (fractional amount of variation of the response data matrix explained by constraints), and applies a correction of the R2 to take into account the number of explanatory variables. This is the statistic that should be reported.
## $r.squared
## [1] 0.3600362
## 
## $adj.r.squared
## [1] 0.3208547

# hypothesis tests #
# testing the first axis (global test)
anova(rda)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance      F Pr(>F)    
## Model     3   6.4807 9.1889  0.001 ***
## Residual 49  11.5193                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda,first=TRUE)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance      F Pr(>F)    
## RDA1      1   5.5244 23.499  0.001 ***
## Residual 49  11.5193                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing all axes sequentially 
anova(rda,by="axis",model="direct",perm.max=9999,step=1000)
## Permutation test for rda under direct model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## RDA1      1   5.5244 23.4991  0.001 ***
## RDA2      1   0.7680  3.2667  0.061 .  
## RDA3      1   0.1883  0.8010  0.512    
## Residual 49  11.5193                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing the individual terms=constraints
anova(rda,by="terms",model="direct",perm.max=9999,step=1000)  # tests terms sequentially, order matters!
## Permutation test for rda under direct model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## logQ      1   1.8477  7.8596  0.001 ***
## logTDN    1   3.9175 16.6639  0.001 ***
## canopy    1   0.7155  3.0434  0.026 *  
## Residual 49  11.5193                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda,by="margin",model="direct",perm.max=9999,step=1000) # tests each term in full model (like drop1() function)
## Permutation test for rda under direct model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## logQ      1   1.0182  4.3312  0.002 ** 
## logTDN    1   2.6904 11.4441  0.001 ***
## canopy    1   0.7155  3.0434  0.026 *  
## Residual 49  11.5193                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

###################
# making triplots #

# again various types of scaling for the plotting step:
# scaling 1 "distance triplot"
# shows similarities between objects/sites in the response matrix because distances among sites reflect their Euclidean distances
# only angles between constraints and responses reflect their correlations (not angles among responses)
# possible conclusions from a "distance triplot":
# Sites that are closer together have more similar communities/water chemistry/DOM composition.
# Species that are closer together occupy more sites in common.

# scaling 2 "correlation triplot"
# shows the effects of constraints/explanatory variables, because all angles between constraints and responses reflect their correlations 
# distances among sites do not reflect their Euclidean distances
# possible conclusions from a "correlation triplot"
# Longer arrows mean this variable strongly drives the variation in the community matrix/water chemistry/DOM composition.
# Arrows pointing in opposite directions have a negative relationship.
# Arrows pointing in the same direction have a positive relationship.

# in both scaling types sites can be projected on constraints and on responses
# factor constraints are shown as centroids instead of arrows, projecting works identical 

# scaling 3 is compromise

# build an RDA scaling type 1 triplot
plot(rda,scaling=1)


(sites<-scores(rda,choices=c(1,2),display="sites",scaling=1)) 
##            RDA1         RDA2
## 1  -0.498460582  0.637194072
## 2  -0.333307201  0.830094727
## 3  -0.235567821  0.585448173
## 4  -0.419169104  0.388569136
## 5  -1.360336512  0.878721129
## 7   0.477322538 -0.075535075
## 8   0.977437916  0.239543485
## 9  -0.267179489  0.018503982
## 10 -0.329982722  0.495341683
## 11  0.249858739  0.269843172
## 12  0.772517930 -0.191121436
## 13 -0.314209914 -0.016205304
## 14  0.814335423 -0.413862255
## 15 -0.572349386 -0.248101430
## 16 -0.631769326  0.146264714
## 17 -0.790465648  0.056555154
## 18 -0.402641922  0.189907030
## 19 -0.311415221 -0.091178758
## 20 -0.081209267 -0.119845614
## 21 -0.439185103 -0.193442177
## 22 -0.249398850 -0.190187243
## 23 -0.400224485 -0.020598570
## 24 -0.717244600 -0.128427989
## 25 -0.142174974 -0.303603256
## 26 -0.003329834 -0.314962678
## 27 -0.092939197 -0.249978470
## 28 -0.118850427 -0.286234700
## 29 -0.084894890 -0.280740961
## 30 -0.337218829 -0.158530543
## 31  0.485881970  0.026354971
## 32  0.290531318  0.070677146
## 33  0.055647198 -0.394945899
## 34  0.062814787  0.276606892
## 35 -0.079826826 -0.188644454
## 36  0.352318114  0.146025624
## 37  0.534072469 -0.103893687
## 38  0.432511681 -0.046075409
## 39  0.555156725  0.066219970
## 40  0.430698723 -0.001927959
## 41  0.761206836 -0.351560915
## 42  0.787558621 -0.228650302
## 43  0.388279575  0.105845304
## 44  0.595488024 -0.222036636
## 45  0.643090631  0.503632310
## 46 -0.447735520 -0.069601236
## 47 -0.109651939  0.026290352
## 48  0.292996889 -0.012505478
## 49 -0.105422063  0.101033121
## 50 -0.245621890 -0.191037745
## 51 -0.246528053 -0.131929977
## 52  0.403636317 -0.289520170
## 53  0.025270757 -0.225210312
## 54 -0.020321586 -0.318575509
## attr(,"const")
## [1] 5.531195

#(lcs<-scores(rda,choices=c(1,2),display="lc",scaling=1)) # fitted/constrained site scores

(species<-scores(rda,choices=c(1,2),display="species",scaling=1)*0.5)
##             RDA1        RDA2
## cond -0.59858199  1.27378225
## pH    0.39342255 -0.19511683
## turb -0.77595131  0.01004908
## Alk  -0.09788046  0.99237637
## TSS  -0.69028937  0.47354883
## PO4  -0.15467198 -0.13922271
## NH4  -0.56926510  0.52603200
## NO2  -0.77353865 -0.43101017
## NO3  -1.00809786 -1.19654781
## Br   -0.86504284  0.28100954
## Cl   -0.75374148 -0.03929891
## Fl   -0.69926370  0.28229280
## SO4  -0.50586602  0.02176366
## Na   -0.67742612  0.17249327
## K    -0.39244249  0.20614834
## Ca   -0.28777201  1.04107246
## Mg   -0.32375066  0.68353998
## TDN  -1.10581086 -1.03980988
## attr(,"const")
## [1] 5.531195

(constraints<-scores(rda,choices=c(1,2),display="bp",scaling=1)*2)
##              RDA1       RDA2
## logQ    0.5732272 -0.2608649
## logTDN -1.0412848 -0.1361117
## canopy  0.5903155 -0.0871586
## attr(,"const")
## [1] 5.531195

perc <- round(100*(summary(rda)$cont$importance[2, 1:2]), 2) #check str(rda) to find such information

plot(sites,asp=1,pch=21,bg=landuse,ylim=c(-1.5,1.5),
     xlab = paste0("RDA1 (", perc[1], "%)"), 
     ylab = paste0("RDA2 (", perc[2], "%)"))
     
Arrows(x0=0,y0=0,x1=constraints[,1],y1=constraints[,2],lwd=1.5,col="blue")
text(constraints[,1:2]*1.1,label=rownames(constraints),pos=4,cex=0.8,col="blue")

Arrows(x0=0,y0=0,x1=species[,1],y1=species[,2],lwd=1,arr.length=0)
text(species[,1:2]*1.1,label=rownames(species),pos=4,cex=0.6)

Unimodal relationships

The classical ordination methods always target dimension reduction

Paliy and Shankar in Molecular Ecology 2016

Correspondence Analysis (CA)

Based on one matrix (usually a species or community matrix).

Considers unimodal responses to (unknown) environmental variables.

Indirect gradient analysis, resulting gradients are synthetic environmental gradients.

PCA on datasets with many zeros produces artefacts (horseshoe effects). Such datasets need ‘unimodal’ methods (e.g. CA) or transform-based PCA or distance-based methods.

CA does not take into account cases where the value of a variable in two different objects is zero.

The basis for CA is weighted averaging from environmental and species tables. If env exists, then this can be done to extract bioindicatory information:

\[ u^*=\frac{y_1x_1+y_2x_2+...+y_nx_n}{y_1+y_2+...+y_n} \] A species optimum \(u^*\) is computed as an abundance-weighted means of a specific environmental variable over all sites at which a specific species is present.

This approach works best when:

Correspondence Analysis (CA)

CA uses a two-way weighted averaging with a theoretical environmental variable iteratively in several steps:

  1. Take arbitrary site scores.
  2. Derive species scores by weighted average of sites scores, for species k (of m): \[ u_k=\sum_{i=1}^n{y_{ki}x_i}/\sum_{i=1}^n{y_{ki}} \]
  3. From the species scores new site scores can be derived, for site i (of n): \[ x_i=\sum_{k=1}^m{y_{ki}u_k}/\sum_{k=1}^n{y_{ki}} \]
  4. Rescaling (standardization) of site and species scores.
  5. Repeat 2-3 several times until stabilisation of site and species scores = first CA axis.
  6. Similar procedure to construct second CA axis (uncorrelated to first).

For studying at home. ;) https://www.davidzeleny.net/anadat-r/doku.php/en:ca_dca

Correspondence Analysis (CA) in R

data(varespec) # a R dataset on vegetation
data(varechem) # soil chemistry
head(varespec) #lichen species
##    Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv Descflex Betupube
## 18     0.55    11.13     0.00     0.00    17.80     0.07     0.00        0
## 15     0.67     0.17     0.00     0.35    12.13     0.12     0.00        0
## 24     0.10     1.55     0.00     0.00    13.47     0.25     0.00        0
## 27     0.00    15.13     2.42     5.92    15.97     0.00     3.70        0
## 23     0.00    12.68     0.00     0.00    23.73     0.03     0.00        0
## 19     0.00     8.92     0.00     2.42    10.28     0.12     0.02        0
##    Vacculig Diphcomp Dicrsp Dicrfusc Dicrpoly Hylosple Pleuschr Polypili
## 18     1.60     2.07   0.00     1.62     0.00      0.0     4.67     0.02
## 15     0.00     0.00   0.33    10.92     0.02      0.0    37.75     0.02
## 24     0.00     0.00  23.43     0.00     1.68      0.0    32.92     0.00
## 27     1.12     0.00   0.00     3.63     0.00      6.7    58.07     0.00
## 23     0.00     0.00   0.00     3.42     0.02      0.0    19.42     0.02
## 19     0.00     0.00   0.00     0.32     0.02      0.0    21.03     0.02
##    Polyjuni Polycomm Pohlnuta Ptilcili Barbhatc Cladarbu Cladrang Cladstel
## 18     0.13     0.00     0.13     0.12     0.00    21.73    21.47     3.50
## 15     0.23     0.00     0.03     0.02     0.00    12.05     8.13     0.18
## 24     0.23     0.00     0.32     0.03     0.00     3.58     5.52     0.07
## 27     0.00     0.13     0.02     0.08     0.08     1.42     7.63     2.55
## 23     2.12     0.00     0.17     1.80     0.02     9.08     9.22     0.05
## 19     1.58     0.18     0.07     0.27     0.02     7.23     4.95    22.08
##    Cladunci Cladcocc Cladcorn Cladgrac Cladfimb Cladcris Cladchlo Cladbotr
## 18     0.30     0.18     0.23     0.25     0.25     0.23     0.00     0.00
## 15     2.65     0.13     0.18     0.23     0.25     1.23     0.00     0.00
## 24     8.93     0.00     0.20     0.48     0.00     0.07     0.10     0.02
## 27     0.15     0.00     0.38     0.12     0.10     0.03     0.00     0.02
## 23     0.73     0.08     1.42     0.50     0.17     1.78     0.05     0.05
## 19     0.25     0.10     0.25     0.18     0.10     0.12     0.05     0.02
##    Cladamau Cladsp Cetreric Cetrisla Flavniva Nepharct Stersp Peltapht Icmaeric
## 18     0.08   0.02     0.02     0.00     0.12     0.02   0.62     0.02        0
## 15     0.00   0.00     0.15     0.03     0.00     0.00   0.85     0.00        0
## 24     0.00   0.00     0.78     0.12     0.00     0.00   0.03     0.00        0
## 27     0.00   0.02     0.00     0.00     0.00     0.00   0.00     0.07        0
## 23     0.00   0.00     0.00     0.00     0.02     0.00   1.58     0.33        0
## 19     0.00   0.00     0.00     0.00     0.02     0.00   0.28     0.00        0
##    Cladcerv Claddefo Cladphyl
## 18        0     0.25        0
## 15        0     1.00        0
## 24        0     0.33        0
## 27        0     0.15        0
## 23        0     1.97        0
## 19        0     0.37        0
head(varechem)
##       N    P     K    Ca    Mg    S    Al   Fe    Mn   Zn  Mo Baresoil Humdepth
## 18 19.8 42.1 139.9 519.4  90.0 32.3  39.0 40.9  58.1  4.5 0.3     43.9      2.2
## 15 13.4 39.1 167.3 356.7  70.7 35.2  88.1 39.0  52.4  5.4 0.3     23.6      2.2
## 24 20.2 67.7 207.1 973.3 209.1 58.1 138.0 35.4  32.1 16.8 0.8     21.2      2.0
## 27 20.6 60.8 233.7 834.0 127.2 40.7  15.4  4.4 132.0 10.7 0.2     18.7      2.9
## 23 23.8 54.5 180.6 777.0 125.8 39.5  24.2  3.0  50.1  6.6 0.3     46.0      3.0
## 19 22.8 40.9 171.4 691.8 151.4 40.8 104.8 17.6  43.6  9.1 0.4     40.5      3.8
##     pH
## 18 2.7
## 15 2.8
## 24 3.0
## 27 2.8
## 23 2.7
## 19 2.7
apply(varespec,1,sum) # approximate 100 (total cover), "absolute" abundance data
##     18     15     24     27     23     19     22     16     28     13     14 
##  89.20  89.82  94.21 125.61  90.46  81.27 109.76  88.52 110.70 101.89  81.65 
##     20     25      7      5      6      3      4      2      9     12     10 
##  64.11  94.06 103.38  94.77 110.90 106.67  84.83 119.13 122.60 119.80 122.37 
##     11     21 
## 112.84  99.17

plot(varechem[,4], varespec[,1])


# correspondence analysis #
# run a CA just based on the species data (unconstrained!)
vare.ca<-cca(X=varespec) # function also used for CCA, but here only one matrix X is supplied

summary(vare.ca,scaling=1)
## 
## Call:
## cca(X = varespec) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total           2.083          1
## Unconstrained   2.083          1
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.5249 0.3568 0.2344 0.19546 0.17762 0.12156 0.11549
## Proportion Explained  0.2520 0.1713 0.1125 0.09383 0.08526 0.05835 0.05544
## Cumulative Proportion 0.2520 0.4233 0.5358 0.62962 0.71489 0.77324 0.82868
##                           CA8     CA9    CA10    CA11    CA12    CA13     CA14
## Eigenvalue            0.08894 0.07318 0.05752 0.04434 0.02546 0.01710 0.014896
## Proportion Explained  0.04269 0.03513 0.02761 0.02129 0.01222 0.00821 0.007151
## Cumulative Proportion 0.87137 0.90650 0.93411 0.95539 0.96762 0.97583 0.982978
##                           CA15     CA16     CA17     CA18     CA19      CA20
## Eigenvalue            0.010160 0.007830 0.006032 0.004008 0.002865 0.0019275
## Proportion Explained  0.004877 0.003759 0.002896 0.001924 0.001375 0.0009253
## Cumulative Proportion 0.987855 0.991614 0.994510 0.996434 0.997809 0.9987341
##                            CA21      CA22      CA23
## Eigenvalue            0.0018074 0.0005864 0.0002434
## Proportion Explained  0.0008676 0.0002815 0.0001168
## Cumulative Proportion 0.9996017 0.9998832 1.0000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##                 CA1       CA2      CA3       CA4        CA5       CA6
## Callvulg  0.0303167 -1.597460  0.11455 -2.894569  0.1376073  2.291129
## Empenigr  0.0751030  0.379305  0.39303  0.023675  0.8568729 -0.400964
## Rhodtome  1.1052309  1.499299  3.04284  0.120106  3.2324306 -0.283510
## Vaccmyrt  1.4614812  1.622935  2.72375  0.231688  0.4604556  0.712538
## Vaccviti  0.1468014  0.313436  0.14696  0.243505  0.6868371 -0.147815
## Pinusylv -0.4820096  0.588517 -0.36020 -0.127094  0.4064754  0.386604
## Descflex  1.5348239  1.218806  1.87562 -0.001340 -1.3136979 -0.070731
## Betupube  0.6694503  1.951826  3.84017  1.389423  7.5959115 -0.244478
## Vacculig -0.0830789 -1.629259  1.05063  0.802648 -0.3058811 -1.625341
## Diphcomp -0.5446464 -1.037570  0.52282  0.940275  0.3682126 -1.082929
## Dicrsp    1.8120408  0.360290 -4.92082  3.088562  1.3867372  0.157815
## Dicrfusc  1.2704743 -0.562978 -0.39718 -2.929542  0.3848272 -2.408710
## Dicrpoly  0.7248118  1.409347  0.80341  1.915549  4.5674148  1.295447
## Hylosple  2.0062408  1.743883  2.27549  0.928884 -3.7648428  2.254851
## Pleuschr  1.3102086  0.583036 -0.01004  0.137298 -1.1216144  0.200422
## Polypili -0.3805097 -1.243904  0.54593  1.477188 -0.7276341 -0.387641
## Polyjuni  1.0133795  0.099043 -2.24697  1.510641  0.7729714 -3.062378
## Polycomm  0.8468241  1.321773  1.13585  1.140723  2.6836594 -0.605038
## Pohlnuta -0.0136453  0.589290 -0.35542  0.135481  0.9369707  0.397246
## Ptilcili  0.4223631  1.598584  3.43474  1.400065  6.3209491  0.198935
## Barbhatc  0.5018348  2.119334  4.57303  1.693188  8.1101807  0.645995
## Cladarbu -0.1531729 -1.483884  0.20024  0.193680  0.0734141  0.358926
## Cladrang -0.5502561 -1.084008  0.40552  0.724060 -0.3357992 -0.335924
## Cladstel -1.4373146  1.077753 -0.44397 -0.375926 -0.2421525  0.004212
## Cladunci  0.8151727 -1.006186 -1.82587 -1.389523  1.6046713  3.675908
## Cladcocc -0.2133215 -0.584429 -0.21434 -0.567886 -0.0003788 -0.145303
## Cladcorn  0.2631227 -0.177858 -0.44464  0.272422  0.3992282 -0.306738
## Cladgrac  0.1956947 -0.311167 -0.23894  0.379013  0.4933026  0.037581
## Cladfimb  0.0009213 -0.161418  0.18463 -0.435908  0.4831233 -0.143751
## Cladcris  0.3373031 -0.470369 -0.05093 -0.823855  0.7182250  0.636140
## Cladchlo -0.6200021  1.207278  0.21889  0.426447  1.9506082  0.120722
## Cladbotr  0.5647242  1.047333  2.65330  0.907734  4.4946805  1.201655
## Cladamau -0.6598144 -1.512880  0.83251  1.577699 -0.0407227 -1.419139
## Cladsp   -0.8209003  0.476164 -0.49752 -0.998241 -0.2393208  0.390785
## Cetreric  0.2458192 -0.689228 -1.68427 -0.131681  0.7439412  2.374535
## Cetrisla -0.3465221  1.362693  0.85897  0.396752  2.7526968  0.396591
## Flavniva -1.4391907 -0.833589 -0.12919  0.007071 -1.4841375  2.956977
## Nepharct  1.6813309  0.199484 -4.33509  2.229917  0.9561223 -5.472858
## Stersp   -0.5172793 -2.280900  0.99775  2.377013 -0.8892757 -1.441228
## Peltapht  0.4035858 -0.043265  0.04538  0.711040  0.1824679 -0.841227
## Icmaeric  0.0378754 -2.419595  0.72135  0.361302 -0.3736424 -2.092136
## Cladcerv -0.9232858 -0.005233 -1.22058  0.305290 -0.8142627  0.414135
## Claddefo  0.5190399 -0.496632 -0.15271 -0.695927  0.9042143  0.909191
## Cladphyl -1.2836161  1.155872 -0.79912 -0.741170 -0.1608002  0.490526
## 
## 
## Site scores (weighted averages of species scores)
## 
##          CA1      CA2       CA3      CA4        CA5      CA6
## 18 -0.108122 -0.53705  0.229574  0.24412  0.1405624 -0.14253
## 15  0.697118 -0.14441 -0.031788 -0.21743 -0.2738522 -0.08146
## 24  0.987603  0.15042 -1.348447  0.80472  0.3095168  0.46773
## 27  0.851765  0.49901  0.443559  0.12277 -0.4814871  0.07589
## 23  0.359881 -0.05608  0.145813  0.15087  0.2405263 -0.17770
## 19  0.003545  0.37017  0.027760  0.06168 -0.1158930 -0.03413
## 22  0.860732 -0.11504  0.110869 -1.02169  0.0772348 -0.60530
## 16  0.636936 -0.33250  0.001120 -0.79797  0.0130769 -0.54049
## 28  1.279352  0.81557  0.670053  0.23137 -0.8929976  0.41783
## 13 -0.195009 -0.80564  0.117686 -0.58286 -0.0007212  0.53071
## 14  0.528532 -0.70420 -0.517771 -0.86836  0.5713441  0.91671
## 20  0.382866 -0.18686 -0.004789  0.10156  0.0458125  0.21087
## 25  0.990715  0.11967 -1.110040  0.44929  0.1885902 -0.70694
## 7  -0.264704 -1.06013  0.334900  0.45973 -0.0326631 -0.19945
## 5  -0.428410 -1.20765  0.374344  0.74970 -0.2596294 -0.30467
## 6  -0.330534 -0.77498  0.130760  0.22391  0.0632686  0.09060
## 3  -0.899601  0.12075 -0.075742  0.03842 -0.1489585 -0.12031
## 4  -0.770294 -0.35351 -0.033779 -0.01795 -0.3007839  0.44303
## 2  -0.992193  0.50319 -0.157505 -0.07070 -0.1065172 -0.09928
## 9  -0.937173  0.78688 -0.258119 -0.19377 -0.0343535 -0.01259
## 12 -0.726413  0.49163 -0.157235 -0.08698 -0.0105774 -0.02801
## 10 -1.002083  0.71239 -0.236526 -0.18643 -0.0231666 -0.04928
## 11 -0.322647 -0.03871 -0.001297  0.09029 -0.1481448  0.06934
## 21  0.259527  0.80746  1.124258  0.36083  1.5437866  0.07051

# As in PCA, the Kaiser-Guttman criterion can be applied to determine the significant axes of a CA.
# Identify the significant axes
ev <- vare.ca$CA$eig
ev[ev > mean(ev)]
##       CA1       CA2       CA3       CA4       CA5       CA6       CA7 
## 0.5249320 0.3567980 0.2344375 0.1954632 0.1776197 0.1215603 0.1154922
n = length(ev)
barplot(ev, main = "Eigenvalues", col = "grey", las = 2)
abline(h = mean(ev), col = "red")
legend("topright", "Average eigenvalue", lwd = 1, col = 2, bty = "n")


# summary(vare.ca,scaling=2)
# again two different types of scaling are possible for biplots

# scaling 1 (distances among sites matter)
# distances among sites approximate their chi^2 distance (=weighted Euclidean distances)
# close sites have similar species abundances
# a site, which is near a specific species, has a high contribution/abundance of that species 

# scaling 2 (relationships among species matter)
# distances among species approximate their chi^2 distance
# close species have similar abundances across sites
# a species, which is near a specific site, is more likely to be found at that site
#layout(matrix(1:2, nrow=1))
plot(vare.ca,scaling=1)

plot(vare.ca,scaling=2)

plot(vare.ca,scaling=3) # a compromise

# for any scaling take care when interpreting species close to origin:
# these are "everywhere" or have optimum right at the origin (i.e., optimum with regard to both axis shown)

# for more controlled plotting
species.scores<-scores(vare.ca,display="species",scaling=2)
site.scores<-scores(vare.ca,display="sites",scaling=2)

plot(site.scores,col="black",pch=21,xlim=c(-2,2),ylim=c(-2,2))
text(species.scores,col="red",label=names(varespec),cex=0.7)
# Species symbols represent their respective optima with regard to the sites, expected abundance of the species decreases with distance from the symbols in all directions (assumption of unimodal distribution). There are no arrows in such a plot because species may have a non-monotonic behaviour across the ordination space!

# post-hoc fitting of an environmental variable
names(varechem)
##  [1] "N"        "P"        "K"        "Ca"       "Mg"       "S"       
##  [7] "Al"       "Fe"       "Mn"       "Zn"       "Mo"       "Baresoil"
## [13] "Humdepth" "pH"
(ef<-envfit(vare.ca,varechem[,12:13],permutations=1999))
## 
## ***VECTORS
## 
##               CA1      CA2     r2 Pr(>r)   
## Baresoil  0.97947 -0.20161 0.2533 0.0565 . 
## Humdepth  0.91602  0.40112 0.4524 0.0040 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1999
plot(ef)

Canonical Correspondence Analysis (CCA)

Unimodal constraint ordination method

Two involved matrices: one dependent, one independent.

In the reciprocal averaging of CA a constraint (e.g. environmental variable) is included:

The result are axes which inform about species-site relationships, but which also have maximized correlation with linear combinations of (environmental) predictors.


Site scores:

Canonical Correspondence Analysis (CCA) in R

# canonical correspondence analysis #
vare.cca<-cca(Y=varespec,X=varechem) 
vare.cca<-cca(varespec~.,varechem) # hypothesis tests need formula interface (don´t ask)

summary(vare.cca,scaling=1)
## 
## Call:
## cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe +      Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.0832      1.000
## Constrained    1.4415      0.692
## Unconstrained  0.6417      0.308
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.16285 0.14213 0.11795 0.08903 0.07029
## Proportion Explained  0.2107 0.1401 0.07817 0.06823 0.05662 0.04274 0.03374
## Cumulative Proportion 0.2107 0.3507 0.42890 0.49713 0.55375 0.59649 0.63023
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.02801 0.01495 0.006382 0.004015 0.003139 0.002955
## Cumulative Proportion 0.65825 0.67319 0.679576 0.683592 0.686730 0.689685
##                          CCA14     CA1     CA2     CA3     CA4     CA5     CA6
## Eigenvalue            0.004733 0.19776 0.14193 0.10117 0.07079 0.05330 0.03330
## Proportion Explained  0.002272 0.09493 0.06813 0.04857 0.03398 0.02559 0.01598
## Cumulative Proportion 0.691958 0.78689 0.85502 0.90359 0.93757 0.96315 0.97914
##                            CA7      CA8      CA9
## Eigenvalue            0.018868 0.015104 0.009488
## Proportion Explained  0.009057 0.007251 0.004554
## Cumulative Proportion 0.988195 0.995446 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3   CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.1628 0.1421 0.11795 0.08903 0.07029
## Proportion Explained  0.3045 0.2024 0.1130 0.0986 0.08183 0.06176 0.04877
## Cumulative Proportion 0.3045 0.5069 0.6198 0.7184 0.80027 0.86203 0.91080
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.04049 0.02160 0.009223 0.005803 0.004536 0.004271
## Cumulative Proportion 0.95128 0.97288 0.982107 0.987910 0.992446 0.996716
##                          CCA14
## Eigenvalue            0.004733
## Proportion Explained  0.003284
## Cumulative Proportion 1.000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##              CCA1     CCA2     CCA3      CCA4     CCA5     CCA6
## Callvulg  0.11374 -1.73246  4.15753  1.844837  3.13741 -1.15626
## Empenigr -0.27373  0.14089  0.09036 -1.134550 -0.40226  0.03525
## Rhodtome -1.59033 -0.11156  0.19187 -2.490433 -0.62292 -1.73616
## Vaccmyrt -1.92827  0.56943  0.75259 -0.244263 -1.65624 -2.05452
## Vaccviti -0.23029  0.22315 -0.13141 -0.960947  0.24441  0.02995
## Pinusylv  0.36674  0.48934  0.55326 -0.726273  0.85051 -0.21227
## Descflex -2.17952  0.50020 -0.40165  1.608949 -1.38617  1.28224
## Betupube -1.07326 -0.41989 -0.20569 -6.388345 -0.62954 -5.60316
## Vacculig  0.77560 -2.19992 -0.93608  0.469588 -2.78966  1.04277
## Diphcomp  0.14991 -1.65300 -1.03898 -1.412059 -0.78833  2.08551
## Dicrsp   -1.28302  0.42862 -4.34136  0.691801  4.43282  1.30777
## Dicrfusc -0.75393 -0.76901  2.04376 -0.684761  0.32655  2.14058
## Dicrpoly -0.79564  0.14903 -2.01239 -3.186680  2.23820 -3.43647
## Hylosple -2.75940  1.46965  0.12345  3.602353 -2.66866 -0.74851
## Pleuschr -1.39625  0.62360 -0.02266  0.817213 -0.19077  0.06281
## Polypili  0.21763 -0.84392 -1.27708 -0.747467 -0.15333  0.16978
## Polyjuni -0.91607  0.38917 -0.87255 -0.891253 -1.78446  1.17847
## Polycomm -1.34974  0.59359 -0.58214 -2.854381 -1.19037 -2.60320
## Pohlnuta -0.01435  0.46778 -0.34834 -0.931562  1.23465 -0.32446
## Ptilcili -0.86964 -0.22650 -0.14520 -5.594842 -0.48392 -5.05263
## Barbhatc -1.04773 -0.42524 -0.29330 -6.830157 -0.50320 -6.88497
## Cladarbu  0.31928 -1.31813 -0.06534  0.138503 -0.11811 -0.26229
## Cladrang  0.57516 -1.14185 -0.60438  0.280955 -0.47617  0.10938
## Cladstel  1.36834  1.29986  0.20555  0.179762 -0.04827  0.09185
## Cladunci -0.34820  0.11797 -0.03422 -1.037583  2.65119 -0.48962
## Cladcocc  0.33121 -0.25213  0.31806 -0.205437  0.09828  0.41903
## Cladcorn -0.34025  0.12974 -0.22432 -0.686052 -0.31883  0.57211
## Cladgrac -0.16429 -0.34432 -0.39566 -0.533216  0.70217 -0.07237
## Cladfimb  0.03022 -0.16993  0.47734 -0.696051 -0.10470 -0.11656
## Cladcris -0.20689  0.02979  1.04812 -1.124294  0.40186 -0.43505
## Cladchlo  0.66964  1.02386 -0.68975 -1.528620  0.49217 -0.75368
## Cladbotr -1.02718 -0.35199  0.48348 -3.528217  0.63524 -4.23041
## Cladamau -0.02415 -2.15363 -1.80591 -1.323302 -1.02050  2.39498
## Cladsp    1.03577  0.72454  0.76099  0.741439  1.75911  0.41843
## Cetreric  0.09754 -0.07199 -1.05941  0.315233  2.75328 -0.58261
## Cetrisla  0.24027  0.64936 -0.12182 -2.346145  0.48511 -2.31098
## Flavniva  1.31684 -1.19677 -1.15320  5.202080  1.07346 -7.81575
## Nepharct -1.15139  0.36798 -1.38414 -0.153781 -3.31081  2.49381
## Stersp    0.18370 -2.37391 -2.38790 -0.009846 -1.07525  1.39790
## Peltapht -0.60047  0.31181  0.12300 -0.899163 -0.76856  0.65021
## Icmaeric  0.26085 -2.83828 -1.06550 -0.409685 -1.20472  1.06913
## Cladcerv  1.06877 -0.10889 -0.78377  3.250753  0.01418 -3.50019
## Claddefo -0.45498 -0.03870  0.60324 -1.497542  0.85219 -0.63272
## Cladphyl  1.51291  2.08492  0.04117 -0.268421  0.27480  0.48796
## 
## 
## Site scores (weighted averages of species scores)
## 
##        CCA1     CCA2      CCA3     CCA4       CCA5     CCA6
## 18  0.11823 -0.57251 -0.164982 -0.22892 -0.1940177  0.07213
## 15 -0.64276 -0.10649  0.169910  0.11432  0.0521025  0.23988
## 24 -0.84786  0.25736 -1.189184  0.14813  1.3580797  0.22853
## 27 -0.99432  0.35227  0.034639  0.28730 -0.4232940 -0.02911
## 23 -0.39622 -0.09941 -0.054725 -0.43892 -0.1038877  0.02098
## 19 -0.07306  0.38585  0.006695 -0.02930 -0.1896126 -0.02464
## 22 -0.72351 -0.26482  0.855780 -0.16216  0.0893297  0.55882
## 16 -0.50071 -0.42518  0.666714 -0.05991  0.1632133  0.51821
## 28 -1.48534  0.62159  0.100450  0.70953 -0.6209907 -0.35786
## 13  0.26728 -0.79352  0.904036  0.45978  0.6372512 -0.27314
## 14 -0.30232 -0.37451  0.439688 -0.39404  0.9278442  0.04663
## 20 -0.36984 -0.13664 -0.135727 -0.13735  0.0942856  0.03259
## 25 -0.85602  0.13551 -0.587776 -0.01017  0.3304823  0.65496
## 7   0.36936 -1.08951 -0.372699  0.05638 -0.4616061  0.05740
## 5   0.44060 -1.21454 -0.658393  0.16630 -0.4226872  0.15976
## 6   0.39218 -0.69770 -0.189710 -0.03141 -0.0990151 -0.05450
## 3   0.88634  0.21282 -0.085773  0.09809 -0.2111364  0.08974
## 4   0.77344 -0.30247 -0.083929  0.80863  0.1228687 -0.94716
## 2   0.93351  0.60859  0.004559  0.01574 -0.1379707  0.08149
## 9   0.86982  0.91296  0.096369 -0.05063  0.0005494  0.01469
## 12  0.67010  0.58561  0.034417 -0.09231 -0.0424654  0.05488
## 10  0.93439  0.83587  0.093851 -0.06296 -0.0540425  0.05003
## 11  0.30814  0.02923 -0.059108  0.09765 -0.0281526 -0.01160
## 21 -0.47641  0.23201  0.003915 -1.44448 -0.2880137 -1.21174
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        CCA1     CCA2      CCA3     CCA4      CCA5     CCA6
## 18 -0.28028 -0.71553 -0.198602 -0.35622 -0.016645  0.28042
## 15 -0.12605  0.26839  0.183427 -0.19962 -0.026309 -0.23569
## 24 -0.57190  0.13619 -1.113919  0.21486  1.130853  0.07844
## 27 -1.12491  0.26289 -0.227400  0.40474 -0.210924  0.14882
## 23 -0.52704  0.05792  0.103917 -0.34088 -0.098759  0.13090
## 19 -0.44851  0.54086  0.013494 -0.37832 -0.048521 -0.27997
## 22 -0.54244 -0.36270  0.612068 -0.02209  0.194629  0.66118
## 16 -0.09855 -0.62779  0.413117 -0.16871 -0.053130 -0.07504
## 28 -1.37258  0.59298  0.200796  0.71143 -0.478758 -0.19023
## 13  0.10953 -0.73196  1.049988  0.47179  0.604495 -0.16296
## 14 -0.09320  0.10867  0.313802 -0.33147  0.232443 -0.11453
## 20 -0.45423  0.04379 -0.082409 -0.42046  0.381970 -0.22780
## 25 -0.59995  0.15944 -0.222688 -0.02782 -0.388701  0.24252
## 7   0.91721 -1.04185 -0.323016  0.13738 -0.567910 -0.03541
## 5   0.06432 -1.09164 -0.636766  0.01507 -0.151542  0.29546
## 6   0.27735 -0.30740 -0.130894  0.02489 -0.019959  0.10057
## 3   0.63365  0.06729 -0.206032  0.05714 -0.365798 -0.04820
## 4   0.56735 -0.42871 -0.189592  0.87651  0.160886 -0.84789
## 2   1.01789  0.50232  0.038999  0.09780 -0.003433  0.21274
## 9   1.01611  0.86649 -0.006136 -0.04395  0.239962  0.19821
## 12  0.29646  0.12958  0.378874 -0.10628  0.044242  0.11423
## 10  0.73605  0.86077 -0.016804  0.04149 -0.158371  0.07950
## 11  0.39119  0.19766 -0.018368 -0.05333 -0.024357 -0.11579
## 21 -0.45499 -0.12585 -0.070007 -1.04926 -0.070611 -0.65096
## 
## 
## Biplot scores for constraining variables
## 
##              CCA1     CCA2      CCA3      CCA4      CCA5       CCA6
## N        -0.14766 -0.28570  0.002715  0.066860 -0.086965  0.0304388
## P        -0.21110  0.31268 -0.065374  0.180759  0.063227 -0.0363527
## K        -0.24254  0.16634  0.145204  0.180743  0.111771 -0.0586722
## Ca       -0.29655  0.22782 -0.015240  0.037048  0.105769  0.0129929
## Mg       -0.28817  0.18393 -0.057371  0.040679  0.170979 -0.0017181
## S        -0.01594  0.22454  0.059879  0.167562  0.205056 -0.0496189
## Al        0.50996 -0.02564  0.015177  0.147400  0.055261 -0.1004202
## Fe        0.43000 -0.04760 -0.016976  0.099139 -0.023974 -0.0332227
## Mn       -0.47852  0.12132  0.045621  0.109904 -0.047628  0.0538484
## Zn       -0.23723  0.18092 -0.112151  0.130337  0.212656 -0.0003564
## Mo        0.13523 -0.05582 -0.063359  0.122240  0.177367 -0.0935487
## Baresoil -0.35558 -0.13762  0.055249 -0.196250  0.057225 -0.1051508
## Humdepth -0.46157  0.10891  0.109612 -0.051173 -0.001117 -0.0153218
## pH        0.32935  0.04056 -0.131693  0.007887 -0.049995 -0.0176316
summary(vare.cca,scaling=2)
## 
## Call:
## cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe +      Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.0832      1.000
## Constrained    1.4415      0.692
## Unconstrained  0.6417      0.308
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.16285 0.14213 0.11795 0.08903 0.07029
## Proportion Explained  0.2107 0.1401 0.07817 0.06823 0.05662 0.04274 0.03374
## Cumulative Proportion 0.2107 0.3507 0.42890 0.49713 0.55375 0.59649 0.63023
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.02801 0.01495 0.006382 0.004015 0.003139 0.002955
## Cumulative Proportion 0.65825 0.67319 0.679576 0.683592 0.686730 0.689685
##                          CCA14     CA1     CA2     CA3     CA4     CA5     CA6
## Eigenvalue            0.004733 0.19776 0.14193 0.10117 0.07079 0.05330 0.03330
## Proportion Explained  0.002272 0.09493 0.06813 0.04857 0.03398 0.02559 0.01598
## Cumulative Proportion 0.691958 0.78689 0.85502 0.90359 0.93757 0.96315 0.97914
##                            CA7      CA8      CA9
## Eigenvalue            0.018868 0.015104 0.009488
## Proportion Explained  0.009057 0.007251 0.004554
## Cumulative Proportion 0.988195 0.995446 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3   CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.1628 0.1421 0.11795 0.08903 0.07029
## Proportion Explained  0.3045 0.2024 0.1130 0.0986 0.08183 0.06176 0.04877
## Cumulative Proportion 0.3045 0.5069 0.6198 0.7184 0.80027 0.86203 0.91080
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.04049 0.02160 0.009223 0.005803 0.004536 0.004271
## Cumulative Proportion 0.95128 0.97288 0.982107 0.987910 0.992446 0.996716
##                          CCA14
## Eigenvalue            0.004733
## Proportion Explained  0.003284
## Cumulative Proportion 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##               CCA1     CCA2      CCA3      CCA4      CCA5      CCA6
## Callvulg  0.075347 -0.93581  1.677742  0.695507  1.077518 -0.345001
## Empenigr -0.181340  0.07610  0.036462 -0.427727 -0.138153  0.010517
## Rhodtome -1.053549 -0.06026  0.077428 -0.938897 -0.213938 -0.518031
## Vaccmyrt -1.277428  0.30759  0.303704 -0.092088 -0.568820 -0.613023
## Vaccviti -0.152563  0.12054 -0.053031 -0.362279  0.083942  0.008938
## Pinusylv  0.242956  0.26432  0.223265 -0.273806  0.292102 -0.063335
## Descflex -1.443872  0.27019 -0.162082  0.606576 -0.476067  0.382590
## Betupube -0.711004 -0.22681 -0.083007 -2.408417 -0.216212 -1.671857
## Vacculig  0.513817 -1.18831 -0.377748  0.177035 -0.958084  0.311138
## Diphcomp  0.099310 -0.89289 -0.419273 -0.532348 -0.270745  0.622270
## Dicrsp   -0.849964  0.23153 -1.751924  0.260810  1.522412  0.390210
## Dicrfusc -0.499460 -0.41539  0.824743 -0.258156  0.112149  0.638702
## Dicrpoly -0.527090  0.08050 -0.812083 -1.201383  0.768689 -1.025365
## Hylosple -1.828026  0.79385  0.049816  1.358093 -0.916528 -0.223338
## Pleuschr -0.924978  0.33684 -0.009146  0.308091 -0.065518  0.018741
## Polypili  0.144172 -0.45586 -0.515356 -0.281796 -0.052660  0.050659
## Polyjuni -0.606869  0.21021 -0.352109 -0.336004 -0.612858  0.351629
## Polycomm -0.894165  0.32063 -0.234919 -1.076106 -0.408823 -0.776736
## Pohlnuta -0.009508  0.25268 -0.140571 -0.351201  0.424031 -0.096811
## Ptilcili -0.576115 -0.12234 -0.058593 -2.109265 -0.166198 -1.507591
## Barbhatc -0.694092 -0.22970 -0.118360 -2.574980 -0.172821 -2.054320
## Cladarbu  0.211517 -0.71201 -0.026366  0.052216 -0.040564 -0.078262
## Cladrang  0.381030 -0.61678 -0.243893  0.105921 -0.163536  0.032637
## Cladstel  0.906486  0.70213  0.082949  0.067771 -0.016579  0.027407
## Cladunci -0.230671  0.06372 -0.013810 -0.391170  0.910527 -0.146092
## Cladcocc  0.219419 -0.13619  0.128350 -0.077450  0.033754  0.125028
## Cladcorn -0.225404  0.07008 -0.090524 -0.258643 -0.109501  0.170706
## Cladgrac -0.108836 -0.18599 -0.159664 -0.201023  0.241156 -0.021594
## Cladfimb  0.020022 -0.09179  0.192626 -0.262413 -0.035959 -0.034780
## Cladcris -0.137056  0.01609  0.422960 -0.423861  0.138016 -0.129810
## Cladchlo  0.443621  0.55305 -0.278345 -0.576292  0.169030 -0.224882
## Cladbotr -0.680481 -0.19013  0.195105 -1.330144  0.218169 -1.262258
## Cladamau -0.015996 -1.16331 -0.728763 -0.498887 -0.350481  0.714608
## Cladsp    0.686166  0.39137  0.307091  0.279524  0.604150  0.124850
## Cetreric  0.064619 -0.03889 -0.427516  0.118844  0.945590 -0.173838
## Cetrisla  0.159171  0.35076 -0.049161 -0.884501  0.166607 -0.689545
## Flavniva  0.872373 -0.64645 -0.465365  1.961193  0.368671 -2.332045
## Nepharct -0.762768  0.19877 -0.558560 -0.057976 -1.137069  0.744096
## Stersp    0.121697 -1.28229 -0.963619 -0.003712 -0.369284  0.417103
## Peltapht -0.397796  0.16843  0.049634 -0.338986 -0.263955  0.194009
## Icmaeric  0.172805 -1.53313 -0.429975 -0.154452 -0.413750  0.319003
## Cladcerv  0.708032 -0.05882 -0.316283  1.225539  0.004871 -1.044377
## Claddefo -0.301412 -0.02090  0.243431 -0.564576  0.292677 -0.188788
## Cladphyl  1.002262  1.12620  0.016613 -0.101195  0.094379  0.145598
## 
## 
## Site scores (weighted averages of species scores)
## 
##       CCA1     CCA2      CCA3     CCA4     CCA5     CCA6
## 18  0.1785 -1.05988 -0.408835 -0.60721 -0.56492  0.24175
## 15 -0.9702 -0.19714  0.421046  0.30324  0.15171  0.80394
## 24 -1.2798  0.47645 -2.946863  0.39292  3.95433  0.76592
## 27 -1.5009  0.65216  0.085837  0.76207 -1.23251 -0.09756
## 23 -0.5981 -0.18404 -0.135611 -1.16425 -0.30249  0.07033
## 19 -0.1103  0.71431  0.016591 -0.07773 -0.55210 -0.08258
## 22 -1.0921 -0.49026  2.120668 -0.43014  0.26010  1.87287
## 16 -0.7558 -0.78712  1.652152 -0.15892  0.47523  1.73677
## 28 -2.2421  1.15075  0.248921  1.88204 -1.80814 -1.19935
## 13  0.4035 -1.46904  2.240249  1.21956  1.85549 -0.91541
## 14 -0.4563 -0.69333  1.089571 -1.04519  2.70161  0.15628
## 20 -0.5583 -0.25296 -0.336340 -0.36433  0.27453  0.10923
## 25 -1.2922  0.25087 -1.456542 -0.02698  0.96227  2.19508
## 7   0.5576 -2.01700 -0.923568  0.14954 -1.34406  0.19237
## 5   0.6651 -2.24847 -1.631533  0.44110 -1.23074  0.53544
## 6   0.5920 -1.29165 -0.470112 -0.08331 -0.28830 -0.18265
## 3   1.3379  0.39399 -0.212551  0.26020 -0.61477  0.30075
## 4   1.1675 -0.55997 -0.207980  2.14490  0.35776 -3.17436
## 2   1.4091  1.12669  0.011297  0.04175 -0.40173  0.27311
## 9   1.3130  1.69016  0.238808 -0.13429  0.00160  0.04923
## 12  1.0115  1.08413  0.085287 -0.24485 -0.12365  0.18392
## 10  1.4105  1.54744  0.232569 -0.16699 -0.15736  0.16768
## 11  0.4651  0.05411 -0.146473  0.25902 -0.08197 -0.03886
## 21 -0.7191  0.42952  0.009702 -3.83149 -0.83861 -4.06109
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        CCA1     CCA2     CCA3     CCA4      CCA5    CCA6
## 18 -0.42308 -1.32466 -0.49215 -0.94489 -0.048464  0.9398
## 15 -0.19026  0.49687  0.45454 -0.52951 -0.076603 -0.7899
## 24 -0.86328  0.25213 -2.76035  0.56993  3.292710  0.2629
## 27 -1.69805  0.48669 -0.56351  1.07358 -0.614147  0.4988
## 23 -0.79557  0.10723  0.25751 -0.90419 -0.287557  0.4387
## 19 -0.67702  1.00130  0.03344 -1.00351 -0.141279 -0.9383
## 22 -0.81881 -0.67147  1.51674 -0.05858  0.566703  2.2159
## 16 -0.14877 -1.16222  1.02373 -0.44751 -0.154699 -0.2515
## 28 -2.07190  1.09778  0.49758  1.88707 -1.394002 -0.6375
## 13  0.16534 -1.35508  2.60193  1.25142  1.760111 -0.5461
## 14 -0.14069  0.20118  0.77762 -0.87922  0.676806 -0.3838
## 20 -0.68566  0.08107 -0.20421 -1.11529  1.112185 -0.7635
## 25 -0.90562  0.29517 -0.55183 -0.07379 -1.131782  0.8128
## 7   1.38453 -1.92877 -0.80045  0.36440 -1.653585 -0.1187
## 5   0.09709 -2.02095 -1.57794  0.03999 -0.441247  0.9902
## 6   0.41866 -0.56908 -0.32436  0.06603 -0.058116  0.3371
## 3   0.95649  0.12458 -0.51056  0.15157 -1.065096 -0.1616
## 4   0.85641 -0.79366 -0.46982  2.32495  0.468453 -2.8417
## 2   1.53650  0.92994  0.09664  0.25941 -0.009995  0.7130
## 9   1.53381  1.60412 -0.01520 -0.11658  0.698700  0.6643
## 12  0.44751  0.23990  0.93887 -0.28191  0.128819  0.3828
## 10  1.11107  1.59354 -0.04164  0.11005 -0.461130  0.2664
## 11  0.59050  0.36592 -0.04552 -0.14145 -0.070919 -0.3881
## 21 -0.68681 -0.23299 -0.17348 -2.78317 -0.205599 -2.1817
## 
## 
## Biplot scores for constraining variables
## 
##              CCA1     CCA2      CCA3     CCA4      CCA5      CCA6
## N        -0.22290 -0.52891  0.006729  0.17735 -0.253216  0.102014
## P        -0.31866  0.57886 -0.162001  0.47947  0.184099 -0.121835
## K        -0.36612  0.30794  0.359824  0.47942  0.325444 -0.196637
## Ca       -0.44764  0.42176 -0.037765  0.09827  0.307969  0.043545
## Mg       -0.43499  0.34051 -0.142169  0.10790  0.497841 -0.005758
## S        -0.02406  0.41570  0.148384  0.44446  0.597063 -0.166296
## Al        0.76978 -0.04747  0.037610  0.39098  0.160905 -0.336554
## Fe        0.64909 -0.08811 -0.042067  0.26297 -0.069806 -0.111345
## Mn       -0.72232  0.22460  0.113052  0.29152 -0.138680  0.180471
## Zn       -0.35810  0.33493 -0.277916  0.34572  0.619191 -0.001195
## Mo        0.20413 -0.10334 -0.157007  0.32424  0.516439 -0.313525
## Baresoil -0.53675 -0.25477  0.136910 -0.52055  0.166621 -0.352409
## Humdepth -0.69673  0.20163  0.271625 -0.13574 -0.003252 -0.051350
## pH        0.49716  0.07509 -0.326341  0.02092 -0.145569 -0.059091
# again different types of scaling are possible for triplots

# hypothesis tests #
# testing the first axis (global test)
anova(vare.cca)
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## Model    14   1.44148 1.4441  0.043 *
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vare.cca,first=TRUE)
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## CCA1      1   0.43887 6.1551  0.025 *
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing all axes sequentially (preceding axes are taken as constraints)
anova(vare.cca,by="axis",model="direct",perm.max=9999,step=1000)
## Permutation test for cca under direct model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## CCA1      1   0.43887 6.1551  0.031 *
## CCA2      1   0.29178 4.0921  0.412  
## CCA3      1   0.16285 2.2839  0.965  
## CCA4      1   0.14213 1.9934  0.979  
## CCA5      1   0.11795 1.6543  0.997  
## CCA6      1   0.08903 1.2486  1.000  
## CCA7      1   0.07029 0.9859  1.000  
## CCA8      1   0.05836 0.8185  1.000  
## CCA9      1   0.03114 0.4367  1.000  
## CCA10     1   0.01329 0.1865  1.000  
## CCA11     1   0.00836 0.1173  1.000  
## CCA12     1   0.00654 0.0917  1.000  
## CCA13     1   0.00616 0.0863  1.000  
## CCA14     1   0.00473 0.0664  1.000  
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing the individual terms=constraints
anova(vare.cca,by="terms",model="direct",perm.max=9999,step=1000)  # tests terms sequentially, order matters!
## Permutation test for cca under direct model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)   
## N         1   0.13001 1.8234  0.057 . 
## P         1   0.17996 2.5240  0.009 **
## K         1   0.13643 1.9134  0.072 . 
## Ca        1   0.08614 1.2081  0.308   
## Mg        1   0.06760 0.9481  0.519   
## S         1   0.17545 2.4606  0.019 * 
## Al        1   0.14797 2.0753  0.035 * 
## Fe        1   0.04982 0.6987  0.690   
## Mn        1   0.05080 0.7124  0.729   
## Zn        1   0.07167 1.0052  0.442   
## Mo        1   0.09585 1.3442  0.235   
## Baresoil  1   0.08683 1.2178  0.261   
## Humdepth  1   0.09502 1.3327  0.223   
## pH        1   0.06794 0.9529  0.469   
## Residual  9   0.64171                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vare.cca,by="margin",model="direct",perm.max=9999,step=1000) # tests each term in full model (like drop1() function)
## Permutation test for cca under direct model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## N         1   0.06302 0.8838  0.555  
## P         1   0.10327 1.4484  0.185  
## K         1   0.10385 1.4565  0.169  
## Ca        1   0.06169 0.8652  0.600  
## Mg        1   0.08830 1.2384  0.263  
## S         1   0.10388 1.4569  0.171  
## Al        1   0.06876 0.9644  0.483  
## Fe        1   0.04765 0.6683  0.772  
## Mn        1   0.08110 1.1375  0.303  
## Zn        1   0.06448 0.9043  0.537  
## Mo        1   0.07761 1.0885  0.377  
## Baresoil  1   0.08966 1.2574  0.257  
## Humdepth  1   0.13548 1.9001  0.047 *
## pH        1   0.06794 0.9529  0.511  
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# quite a lot of variables in the constraining matrix, maybe selection would be adequate
# --> function ordistep()

###################
# making triplots #

# again various types of scaling for the plotting step:
# scaling 1 "distance triplot"
# sites can be projected on constraints
# sites close to centroid of factor constraint are more likely to possess the specific state (factor level)
# distances among sites reflect their Chi^2 distances

# scaling 2 "correlation triplot"
# species can be projected on constraints (to give their optimum)
# species close to centroid of factor constraint are more likely to be found in the respective sites
# distances among sites do not reflect their Chi^2 distances

# scaling 3 is compromise

plot(vare.cca,scaling=1,display=c("species","sites"))
plot(vare.cca,scaling=1,display=c("species","sites","bp"))

plot(vare.cca,scaling=2)

plot(vare.cca,scaling=3) # a compromise

# for any scaling take care when interpreting species close to origin:
# these are "everywhere" or have optimum right at the origin (i.e., optimum with regard to both axis shown)

# for more controlled plotting compute scores individually
(species<-scores(vare.cca,display="species",scaling=2))
##                  CCA1        CCA2
## Callvulg  0.075346828 -0.93580864
## Empenigr -0.181339634  0.07610159
## Rhodtome -1.053549303 -0.06026078
## Vaccmyrt -1.277428382  0.30758571
## Vaccviti -0.152563159  0.12053851
## Pinusylv  0.242955730  0.26432438
## Descflex -1.443872268  0.27018775
## Betupube -0.711004196 -0.22680891
## Vacculig  0.513817276 -1.18831302
## Diphcomp  0.099309627 -0.89289046
## Dicrsp   -0.849963650  0.23152562
## Dicrfusc -0.499459808 -0.41538963
## Dicrpoly -0.527089621  0.08050066
## Hylosple -1.828025570  0.79385044
## Pleuschr -0.924977860  0.33684271
## Polypili  0.144171903 -0.45585592
## Polyjuni -0.606869478  0.21021273
## Polycomm -0.894165298  0.32063347
## Pohlnuta -0.009508326  0.25267888
## Ptilcili -0.576115080 -0.12234488
## Barbhatc -0.694091803 -0.22969682
## Cladarbu  0.211516504 -0.71200604
## Cladrang  0.381029817 -0.61678185
## Cladstel  0.906485693  0.70213449
## Cladunci -0.230670950  0.06372124
## Cladcocc  0.219418695 -0.13619040
## Cladcorn -0.225403730  0.07008025
## Cladgrac -0.108836318 -0.18598627
## Cladfimb  0.020022402 -0.09178764
## Cladcris -0.137056217  0.01608948
## Cladchlo  0.443621232  0.55305019
## Cladbotr -0.680480562 -0.19013376
## Cladamau -0.015995645 -1.16331056
## Cladsp    0.686166357  0.39136838
## Cetreric  0.064619228 -0.03888872
## Cetrisla  0.159171117  0.35076227
## Flavniva  0.872372656 -0.64644751
## Nepharct -0.762767923  0.19876924
## Stersp    0.121696721 -1.28229477
## Peltapht -0.397796035  0.16842583
## Icmaeric  0.172805490 -1.53313439
## Cladcerv  0.708032101 -0.05881737
## Claddefo -0.301411810 -0.02090192
## Cladphyl  1.002261821  1.12619548

(lcs<-scores(vare.cca,display="lc",scaling=2)) # fitted site scores
##           CCA1        CCA2
## 18 -0.42308410 -1.32465581
## 15 -0.19026462  0.49687345
## 24 -0.86328447  0.25212611
## 27 -1.69805066  0.48669129
## 23 -0.79556975  0.10723446
## 19 -0.67702233  1.00129602
## 22 -0.81881442 -0.67146524
## 16 -0.14876518 -1.16222054
## 28 -2.07190438  1.09778134
## 13  0.16533766 -1.35507570
## 14 -0.14068713  0.20117513
## 20 -0.68566132  0.08107181
## 25 -0.90561895  0.29516717
## 7   1.38453057 -1.92876582
## 5   0.09708596 -2.02094655
## 6   0.41865544 -0.56908384
## 3   0.95649196  0.12457945
## 4   0.85641341 -0.79366326
## 2   1.53649728  0.92993610
## 9   1.53380971  1.60412304
## 12  0.44750585  0.23989837
## 10  1.11106621  1.59354430
## 11  0.59049933  0.36592318
## 21 -0.68681020 -0.23299449

(sites<-scores(vare.cca,display="sites",scaling=2)) # unfitted site scores
##          CCA1        CCA2
## 18  0.1784733 -1.05988423
## 15 -0.9702382 -0.19713866
## 24 -1.2798478  0.47644978
## 27 -1.5009195  0.65215591
## 23 -0.5980933 -0.18403623
## 19 -0.1102881  0.71431421
## 22 -1.0921288 -0.49026245
## 16 -0.7558244 -0.78712482
## 28 -2.2421137  1.15075172
## 13  0.4034539 -1.46904155
## 14 -0.4563466 -0.69332916
## 20 -0.5582740 -0.25295922
## 25 -1.2921584  0.25086578
## 7   0.5575516 -2.01700439
## 5   0.6650822 -2.24846553
## 6   0.5919915 -1.29165187
## 3   1.3379214  0.39399435
## 4   1.1674982 -0.55996655
## 2   1.4091299  1.12668910
## 9   1.3129824  1.69015916
## 12  1.0115068  1.08413112
## 10  1.4104517  1.54743675
## 11  0.4651334  0.05410696
## 21 -0.7191331  0.42951771

(constraints<-scores(vare.cca,choices=c(1,2),display="bp",scaling=2))
##                 CCA1        CCA2
## N        -0.22289549 -0.52891359
## P        -0.31865680  0.57886019
## K        -0.36611630  0.30794063
## Ca       -0.44763820  0.42176452
## Mg       -0.43499067  0.34051494
## S        -0.02405594  0.41569739
## Al        0.76977830 -0.04747443
## Fe        0.64908911 -0.08811470
## Mn       -0.72232495  0.22459941
## Zn       -0.35809796  0.33492792
## Mo        0.20413234 -0.10334347
## Baresoil -0.53674724 -0.25477013
## Humdepth -0.69673051  0.20163085
## pH        0.49715734  0.07508792

plot(sites,col="black",pch=21,xlim=c(-2,2),ylim=c(-2,2))
#text(species,col="red",label=names(varespec),cex=0.7)
text(lcs,col="red",label=names(varespec),cex=0.7)
Arrows(x0=0,y0=0,x1=constraints[,1],y1=constraints[,2],lwd=1.5,col="blue")
text(constraints[,1:2]*1.1,label=rownames(constraints),pos=4,cex=0.8,col="blue")