Carga de paquetes
library(ade4)
library(adegraphics)
##
## Attaching package: 'adegraphics'
## The following objects are masked from 'package:ade4':
##
## kplotsepan.coa, s.arrow, s.class, s.corcircle, s.distri, s.image, s.label, s.logo, s.match,
## s.traject, s.value, table.value, triangle.class
library(adespatial)
##
## Attaching package: 'adespatial'
## The following object is masked from 'package:ade4':
##
## multispati
library(cocorresp)
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-4
##
## Attaching package: 'cocorresp'
## The following object is masked from 'package:ade4':
##
## coinertia
library(vegan)
library(vegan3d)
library(ape)
##
## Attaching package: 'ape'
## The following object is masked from 'package:adegraphics':
##
## zoom
library(MASS)
library(ellipse)
##
## Attaching package: 'ellipse'
## The following object is masked from 'package:graphics':
##
## pairs
library(FactoMineR)
##
## Attaching package: 'FactoMineR'
## The following object is masked from 'package:ade4':
##
## reconst
library(rrcov)
## Loading required package: robustbase
## Scalable Robust Estimators with High Breakdown Point (version 1.4-7)
library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.0 ✔ purrr 0.3.2
## ✔ tibble 2.1.1 ✔ dplyr 0.8.0.1
## ✔ tidyr 0.8.3 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
Funciones epecesciales
source("hcoplot.R")
source("triplot.rda.R")
source("plot.lda.R")
source("polyvars.R")
source("screestick.R")
Lectura de las bases de datos.
# Lectura de las bases de datos.
ambient <- read.csv2("ambientales.csv", enc="latin1", row.names=1)
peces <- read.csv2("peces.csv", enc="latin1", row.names=1)
locs <- read.csv2("localidades.csv", enc="latin1", row.names=1)
nomEpecesc <- read.csv2("nombresdeespecies.csv", enc="latin1")
nomVarAmbient <- read.csv2("Nombresdevariablesambientales.csv", enc="latin1")
Adecuación de la base de datos.
# Adecuación
# Quitar el sitio que no tiene ninguna epecescie
sumEsp <- apply(peces, 1, sum)
peces <- peces[sumEsp != 0, ]
ambient <- ambient[sumEsp != 0, ]
locs <- locs[sumEsp != 0, ]
Análisis de correspondencia canónica.
(peces.cca <- cca(peces ~ ., ambient))
## Call: cca(formula = peces ~ das + alt + slo + flo + pH + har + pho + nit + amm + oxy + bdo, data =
## ambient)
##
## Inertia Proportion Rank
## Total 1.1669 1.0000
## Constrained 0.8377 0.7179 11
## Unconstrained 0.3292 0.2821 17
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10 CCA11
## 0.5345 0.1227 0.0688 0.0489 0.0275 0.0130 0.0098 0.0054 0.0035 0.0022 0.0016
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.10182 0.05329 0.05021 0.03452 0.03000 0.01456 0.01226 0.00887
## (Showing 8 of 17 unconstrained eigenvalues)
summary(peces.cca)
##
## Call:
## cca(formula = peces ~ das + alt + slo + flo + pH + har + pho + nit + amm + oxy + bdo, data = ambient)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 1.1669 1.0000
## Constrained 0.8377 0.7179
## Unconstrained 0.3292 0.2821
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## Eigenvalue 0.5345 0.1227 0.06876 0.04886 0.02746 0.01295 0.009775 0.00545 0.003523 0.002170
## Proportion Explained 0.4580 0.1052 0.05892 0.04187 0.02353 0.01110 0.008377 0.00467 0.003019 0.001859
## Cumulative Proportion 0.4580 0.5632 0.62211 0.66399 0.68752 0.69862 0.706997 0.71167 0.714686 0.716545
## CCA11 CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9
## Eigenvalue 0.001596 0.10182 0.05329 0.05021 0.03452 0.03000 0.01456 0.01226 0.008865 0.007317
## Proportion Explained 0.001368 0.08726 0.04567 0.04303 0.02958 0.02571 0.01248 0.01050 0.007597 0.006270
## Cumulative Proportion 0.717914 0.80517 0.85084 0.89387 0.92345 0.94917 0.96165 0.97215 0.979744 0.986015
## CA10 CA11 CA12 CA13 CA14 CA15 CA16 CA17
## Eigenvalue 0.004570 0.003873 0.002968 0.002148 0.001493 0.0008308 0.0003886 4.723e-05
## Proportion Explained 0.003917 0.003319 0.002544 0.001841 0.001280 0.0007120 0.0003330 4.048e-05
## Cumulative Proportion 0.989931 0.993251 0.995794 0.997635 0.998915 0.9996265 0.9999595 1.000e+00
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10
## Eigenvalue 0.5345 0.1227 0.06876 0.04886 0.02746 0.01295 0.009775 0.005450 0.003523 0.00217
## Proportion Explained 0.6380 0.1465 0.08207 0.05833 0.03278 0.01546 0.011668 0.006505 0.004205 0.00259
## Cumulative Proportion 0.6380 0.7845 0.86656 0.92489 0.95767 0.97313 0.984794 0.991299 0.995504 0.99809
## CCA11
## Eigenvalue 0.001596
## Proportion Explained 0.001906
## 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
## Cogo -1.276782 1.35155 -0.012375 0.22531 -0.179903 0.094379
## Satr -1.525689 -0.43466 -0.328892 0.23630 0.144768 -0.007892
## Phph -1.221136 -0.19104 0.007043 -0.05474 0.018108 -0.045958
## Neba -0.966876 -0.29832 0.078083 -0.12854 0.052918 0.061401
## Thth -1.324767 1.38168 -0.151752 0.48930 -0.180045 -0.293303
## Teso -0.895169 1.18994 0.050588 -0.06064 -0.162988 0.235364
## Chna 0.494080 0.18453 0.118825 -0.10125 0.356934 0.081358
## Chto 0.151333 0.43227 0.125378 -0.53252 0.352606 -0.021676
## Lele -0.119548 0.03933 0.260802 -0.15326 -0.191334 0.051822
## Lece -0.003793 -0.08117 0.386934 0.01045 -0.255979 0.006033
## Baba 0.305425 0.28556 -0.108575 -0.18426 0.185070 0.114957
## Spbi 0.354748 0.31834 -0.044857 -0.41401 0.222240 -0.351709
## Gogo 0.262245 0.02464 0.094701 0.07675 0.113405 0.065630
## Eslu 0.166490 -0.28428 0.007552 0.02723 -0.130462 -0.084917
## Pefl 0.144865 -0.16017 0.054768 -0.26509 -0.145556 -0.168154
## Rham 0.578999 0.07021 -0.304787 -0.12841 0.100704 -0.052385
## Legi 0.620722 0.04340 -0.207869 -0.01401 0.030758 0.011737
## Scer 0.523613 -0.14218 -0.008426 0.13359 -0.163790 -0.216274
## Cyca 0.592932 0.06934 -0.382972 -0.05516 -0.071659 0.034195
## Titi 0.295852 -0.18634 -0.010118 -0.14245 -0.070774 0.152197
## Abbr 0.700227 -0.04107 -0.444469 0.08575 -0.065584 0.126693
## Icme 0.779693 -0.08488 -0.700727 0.17060 -0.349148 -0.167956
## Acce 0.758681 -0.07352 -0.075256 0.26356 -0.009589 0.100114
## Ruru 0.380007 -0.14619 0.392859 -0.05620 -0.128640 -0.006855
## Blbj 0.744602 -0.05553 -0.330245 0.11600 -0.050556 0.140966
## Alal 0.670437 0.01773 0.410319 0.53876 0.271721 -0.064363
## Anan 0.651162 0.01894 -0.384970 -0.05095 -0.032570 -0.058778
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## 1 -2.85448 -3.542243 -4.78349 4.83584 5.271517 -0.6094
## 2 -2.40317 -2.602684 -1.67506 0.98386 2.897990 -0.2515
## 3 -2.15182 -2.498007 -1.10107 0.37389 2.158658 -0.2276
## 4 -1.41353 -2.063451 -0.20347 -0.48937 0.292725 -0.8701
## 5 -0.20080 -1.272187 1.69105 -1.16361 -3.297892 -1.4784
## 6 -1.14732 -1.793445 0.68997 -0.76727 -0.632502 1.0056
## 7 -2.04292 -2.277225 -0.52553 0.22641 1.396340 0.6832
## 9 -0.31069 -1.317289 3.88257 -1.10418 -4.391402 1.6171
## 10 -1.37143 -1.426835 1.45628 -1.03158 -0.916454 1.2974
## 11 -2.21687 0.211582 -0.86950 2.02128 -0.186283 -2.3797
## 12 -2.24382 0.567984 -1.00620 2.14659 0.064873 -1.5851
## 13 -2.36193 2.360284 -1.40230 2.63669 -0.585878 -1.4910
## 14 -1.91679 2.553732 -0.73172 2.02824 -1.054392 -0.2921
## 15 -1.40724 1.948920 0.38348 0.44007 -1.601560 3.3146
## 16 -0.76771 1.491747 0.73233 -2.02643 -0.161016 3.1721
## 17 -0.32180 1.012062 0.60116 -2.04569 1.611913 -1.7428
## 18 -0.07148 0.812691 0.63263 -1.62433 1.139577 -1.7708
## 19 0.19279 0.099176 0.90109 -1.25404 1.701956 0.4343
## 20 0.59825 0.003113 0.64635 -0.59114 1.200796 -0.4261
## 21 0.70530 -0.074133 -0.03580 -0.19148 0.618312 0.1850
## 22 0.76014 -0.073194 -0.23365 -0.06468 -0.005157 0.9481
## 23 0.80314 -0.390997 5.81927 5.27872 1.445825 -2.5009
## 24 0.96586 -0.171744 3.04080 4.37524 2.477896 0.9084
## 25 0.72138 -0.433817 3.15322 3.71187 0.536863 -1.4876
## 26 0.82289 -0.272389 0.06034 1.06446 0.175249 0.7865
## 27 0.81044 -0.213067 -0.44897 0.48279 -0.250576 1.4377
## 28 0.84855 -0.199638 -0.85082 0.48169 -0.742223 0.6115
## 29 0.65716 0.130297 -0.84172 0.01666 -0.336119 -0.6828
## 30 0.86319 -0.098907 -1.18924 -0.16484 -0.840025 -1.8059
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## 1 -3.32694 -2.89515 -3.29235 2.20626 2.3060 1.19734
## 2 -1.14185 -3.48021 -0.43333 -1.28383 2.3367 -0.31102
## 3 -1.93997 -2.25021 -0.99858 0.51805 0.6466 0.28492
## 4 -1.52691 -1.99021 -0.02728 -0.74518 -0.0490 -0.11014
## 5 -0.88266 -1.27553 2.02571 -0.61079 -2.4324 -0.61528
## 6 -1.64941 -2.17812 -0.45378 0.86544 1.3469 -0.17736
## 7 -1.99263 -1.43304 -0.43014 0.89711 -1.5779 0.09591
## 9 0.13595 -1.35685 2.65394 -2.09948 -2.8253 -0.16589
## 10 -1.30135 -0.70616 0.02396 -0.02439 -0.2340 0.56841
## 11 -2.14664 1.68069 -0.44700 -0.32485 -1.5256 0.75702
## 12 -1.31891 0.71524 0.03340 0.39201 1.2446 -1.01908
## 13 -1.57343 1.39536 -0.04734 0.87388 -0.3098 -0.48242
## 14 -1.96408 1.99860 -0.50752 1.54706 -0.5462 -0.09997
## 15 -1.38328 1.80655 -0.17564 0.15301 -0.7491 0.98224
## 16 -0.44297 0.57049 0.29502 -0.89002 -0.1082 1.20984
## 17 -0.16066 1.13332 0.97402 -0.97270 0.5673 -1.34255
## 18 -0.13172 0.68153 0.77466 -1.17474 0.4923 -0.06227
## 19 -0.31662 0.78989 0.20453 -0.49484 1.1835 0.33682
## 20 0.27681 0.25870 0.29757 -0.98233 1.5082 -0.61280
## 21 0.72426 -0.11812 0.28673 -0.81779 0.9103 0.70093
## 22 0.42432 0.18792 -0.10587 0.69590 0.2200 1.32655
## 23 0.08867 0.16111 2.12318 6.90639 1.5586 -1.01857
## 24 1.67581 0.01813 3.09214 1.76719 -0.4168 -0.47757
## 25 0.82273 -0.77842 3.79368 4.54032 1.8147 -2.08635
## 26 0.91779 -0.27125 0.64570 0.49523 -0.2030 1.35326
## 27 0.99628 -0.45687 -0.45261 0.38428 -0.1698 1.18499
## 28 0.96393 -0.23682 -0.80926 0.09871 -0.3918 0.49065
## 29 0.54853 0.00564 -1.12139 -0.08810 -0.1117 -1.83093
## 30 0.89306 0.01746 -1.10304 0.25786 -1.0940 -0.84974
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## das 0.8583 0.22288 -0.38283 0.181558 -0.06361 -0.118244
## alt -0.7941 -0.52544 0.22871 -0.119503 -0.11196 0.003147
## slo -0.6912 -0.39122 0.08569 -0.001661 0.07428 0.269199
## flo 0.6870 0.21431 -0.49082 0.154243 -0.15717 -0.376725
## pH -0.1547 0.26603 -0.18601 0.170969 -0.34080 0.310399
## har 0.5447 0.48531 -0.18551 0.172547 -0.43485 -0.379671
## pho 0.4259 -0.03050 0.38493 0.589627 0.02831 -0.110972
## nit 0.6839 0.20008 0.16675 0.143281 0.25238 0.115377
## amm 0.4114 -0.09372 0.53787 0.473945 0.21177 -0.174246
## oxy -0.7796 0.35802 -0.26759 -0.280475 0.23527 -0.052069
## bdo 0.4550 -0.17355 0.50624 0.633983 -0.05423 -0.025172
RsquareAdj(peces.cca)
## $r.squared
## [1] 0.7179135
##
## $adj.r.squared
## [1] 0.5495456
- sp: epecescies.
- sites: sitios.
- lc: linear onstraints
- bp: flechas biplot
- cn: centroides de variables explicativas categóricas.
- reg: coeficientes de regresion.
op <- par(no.readonly = TRUE)
par(mfrow = c(2, 1))
plot(peces.cca,
scaling = 1,
display = c("sp", "lc", "cn"),
main = "Triplot CCA peces ~ ambient - escalamiento 1"
)
text(-4, 3.9, "a", cex = 1.5)
plot(peces.cca,
scaling = 2,
display = c("sp", "lc", "cn"),
main = "Triplot CCA peces ~ ambient - escalamiento 2")
text(-4, 1.6, "b", cex = 1.5)
par(op)
plot(peces.cca,
scaling = 1,
display = c("lc", "cn"),
main = "Triplot CCA peces ~ ambient - escalamiento 1"
)
plot(peces.cca,
scaling = 2,
display = c("sp", "cn"),
main = "Triplot CCA peces ~ ambient - escalamiento 2")
anova(peces.cca, permutations = how(nperm = 999))
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = peces ~ das + alt + slo + flo + pH + har + pho + nit + amm + oxy + bdo, data = ambient)
## Df ChiSquare F Pr(>F)
## Model 11 0.83774 3.9332 0.001 ***
## Residual 17 0.32917
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(peces.cca, by = "axis", permutations = how(nperm = 999))
## Permutation test for cca under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = peces ~ das + alt + slo + flo + pH + har + pho + nit + amm + oxy + bdo, data = ambient)
## Df ChiSquare F Pr(>F)
## CCA1 1 0.53449 27.6037 0.001 ***
## CCA2 1 0.12271 6.3372 0.001 ***
## CCA3 1 0.06876 3.5509 0.256
## CCA4 1 0.04886 2.5236 0.631
## CCA5 1 0.02746 1.4183 0.981
## CCA6 1 0.01295 0.6688 1.000
## CCA7 1 0.00977 0.5048 1.000
## CCA8 1 0.00545 0.2815 1.000
## CCA9 1 0.00352 0.1819 1.000
## CCA10 1 0.00217 0.1121 1.000
## CCA11 1 0.00160 0.0824 1.000
## Residual 17 0.32917
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cca.step.forward <-
ordistep(cca(peces ~ 1, data = ambient),
scope = formula(peces.cca),
direction = "forward",
permutations = how(nperm = 199))
##
## Start: peces ~ 1
##
## Df AIC F Pr(>F)
## + das 1 98.640 14.7293 0.005 **
## + alt 1 99.999 12.8184 0.005 **
## + oxy 1 100.885 11.6211 0.005 **
## + flo 1 103.362 8.4593 0.005 **
## + slo 1 103.406 8.4052 0.005 **
## + nit 1 103.912 7.7930 0.005 **
## + har 1 105.857 5.5355 0.005 **
## + amm 1 107.962 3.2577 0.005 **
## + pho 1 107.941 3.2800 0.010 **
## + bdo 1 107.219 4.0431 0.015 *
## + pH 1 110.469 0.7517 0.505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: peces ~ das
##
## Df AIC F Pr(>F)
## + oxy 1 94.773 5.8296 0.005 **
## + alt 1 96.941 3.5371 0.005 **
## + flo 1 97.670 2.8043 0.015 *
## + bdo 1 97.455 3.0181 0.020 *
## + amm 1 97.910 2.5670 0.020 *
## + har 1 98.161 2.3202 0.025 *
## + nit 1 98.643 1.8532 0.085 .
## + pho 1 98.732 1.7680 0.100 .
## + pH 1 99.489 1.0527 0.300
## + slo 1 99.496 1.0457 0.365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: peces ~ das + oxy
##
## Df AIC F Pr(>F)
## + alt 1 92.191 4.2798 0.005 **
## + bdo 1 93.471 3.0148 0.010 **
## + pho 1 94.520 2.0195 0.050 *
## + har 1 94.534 2.0071 0.055 .
## + flo 1 94.775 1.7834 0.090 .
## + amm 1 95.145 1.4434 0.170
## + nit 1 95.204 1.3903 0.180
## + slo 1 95.350 1.2577 0.250
## + pH 1 95.971 0.7012 0.730
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: peces ~ das + oxy + alt
##
## Df AIC F Pr(>F)
## + bdo 1 90.110 3.6264 0.005 **
## + pho 1 91.527 2.3088 0.020 *
## + har 1 91.608 2.2354 0.035 *
## + flo 1 92.461 1.4747 0.110
## + amm 1 92.805 1.1742 0.310
## + slo 1 93.249 0.7917 0.560
## + pH 1 93.231 0.8075 0.580
## + nit 1 93.725 0.3885 0.935
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: peces ~ das + oxy + alt + bdo
##
## Df AIC F Pr(>F)
## + har 1 89.567 2.1075 0.030 *
## + flo 1 90.400 1.3972 0.180
## + pH 1 90.913 0.9690 0.435
## + slo 1 91.036 0.8678 0.550
## + amm 1 91.326 0.6298 0.715
## + pho 1 91.476 0.5086 0.875
## + nit 1 91.581 0.4236 0.930
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: peces ~ das + oxy + alt + bdo + har
##
## Df AIC F Pr(>F)
## + pH 1 90.279 0.9997 0.375
## + flo 1 90.606 0.7417 0.605
## + slo 1 90.683 0.6811 0.730
## + amm 1 90.764 0.6180 0.750
## + pho 1 90.887 0.5220 0.830
## + nit 1 91.212 0.2714 0.980
peces.cca.pars <- cca(peces ~ das + oxy + alt + bdo + har, data = ambient)
anova(peces.cca.pars, permutations = how(nperm = 999))
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = peces ~ das + oxy + alt + bdo + har, data = ambient)
## Df ChiSquare F Pr(>F)
## Model 5 0.74784 8.2088 0.001 ***
## Residual 23 0.41907
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(peces.cca.pars, permutations = how(nperm = 999), by = "axis")
## Permutation test for cca under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = peces ~ das + oxy + alt + bdo + har, data = ambient)
## Df ChiSquare F Pr(>F)
## CCA1 1 0.50990 27.9849 0.001 ***
## CCA2 1 0.11315 6.2101 0.001 ***
## CCA3 1 0.06272 3.4424 0.008 **
## CCA4 1 0.04036 2.2150 0.061 .
## CCA5 1 0.02172 1.1918 0.269
## Residual 23 0.41907
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(peces.cca.pars)
## $r.squared
## [1] 0.640873
##
## $adj.r.squared
## [1] 0.5661042
vif.cca(peces.cca)
## das alt slo flo pH har pho nit amm oxy
## 168.535026 37.145185 5.081439 74.270388 2.078310 4.410247 9.378410 11.329111 11.723807 15.229089
## bdo
## 9.913511
vif.cca(peces.cca.pars)
## das oxy alt bdo har
## 14.846811 6.736614 8.644642 4.194945 2.686473
ordirgl(peces.cca.pars,
type = "t",
scaling = 1)
Análisis de correlación canónica.
ambientQuim <- ambient %>%
mutate(logpho = log(pho),
rqnit = sqrt(nit),
logamm = log1p(amm),
logbdo = log(bdo)) %>%
select(logpho, rqnit, logamm, logbdo, har, oxy)
ambientTopo <- ambient %>%
mutate(logElev = log(alt),
logSlo = log(slo),
rqDas = log(das)) %>%
select(logElev, logSlo, rqDas)
cancor(ambientQuim, ambientTopo)
## $cor
## [1] 0.9293904 0.6362018 0.5204193
##
## $xcoef
## [,1] [,2] [,3] [,4] [,5] [,6]
## logpho -0.015685435 0.229875370 0.054931021 3.938770e-03 0.2106656384 -6.502967e-02
## rqnit -0.016343959 -0.062293173 -0.008680533 -2.216213e-02 0.0309164408 2.061974e-02
## logamm -0.003401055 0.079285337 0.137086506 -1.375295e-02 -0.1612218076 -1.884567e-01
## logbdo 0.024926300 -0.159825578 -0.235988986 4.086169e-01 0.0004811065 -3.905072e-01
## har -0.006390052 0.003349485 -0.006171565 -5.035175e-04 -0.0151480300 9.204106e-05
## oxy -0.000480441 0.003311146 -0.002039773 1.092565e-06 0.0008506728 -2.331344e-02
##
## $ycoef
## [,1] [,2] [,3]
## logElev 0.13499570 0.6582992 -0.1542593
## logSlo 0.02135998 -0.5559215 -0.5893563
## rqDas -0.06661458 0.1092621 -0.1628288
##
## $xcenter
## logpho rqnit logamm logbdo har oxy
## 3.213364 11.822361 2.002743 3.713075 85.827586 94.724138
##
## $ycenter
## logElev logSlo rqDas
## 5.9968391 0.9349566 4.6187570
quim.topo.ccora <-
CCorA(ambientQuim, ambientTopo,
stand.Y = TRUE,
stand.X = TRUE,
permutations = how(nperm = 999))
quim.topo.ccora
##
## Canonical Correlation Analysis
##
## Call:
## CCorA(Y = ambientQuim, X = ambientTopo, stand.Y = TRUE, stand.X = TRUE, permutations = how(nperm = 999))
##
## Y X
## Matrix Ranks 6 3
##
## Pillai's trace: 1.539355
##
## Significance of Pillai's trace:
## from F-distribution: 2.8305e-05
## based on permutations: 0.001
## Permutation: free
## Number of permutations: 999
##
## CanAxis1 CanAxis2 CanAxis3
## Canonical Correlations 0.92939 0.63620 0.5204
##
## Y | X X | Y
## RDA R squares 0.51817 0.7680
## adj. RDA R squares 0.46035 0.7047
biplot(quim.topo.ccora, plot.type = "biplot")
biplot(quim.topo.ccora)
Análisis discriminante.
peces.hel <- decostand(peces, "hellinger")
gr <- cutree(hclust(vegdist(peces.hel, "euclidean"), "ward.D2"), k = 4)
ambient3 <- ambient %>%
select(alt, oxy, bdo)
ambient3.d1 <- dist(ambient3)
(ambient.MHV <- betadisper(ambient3.d1, gr))
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = ambient3.d1, group = gr)
##
## No. of Positive Eigenvalues: 3
## No. of Negative Eigenvalues: 0
##
## Average distance to median:
## 1 2 3 4
## 198.76 217.31 40.49 21.87
##
## Eigenvalues for PCoA axes:
## PCoA1 PCoA2 PCoA3
## 2044986 44023 3149
anova(ambient.MHV)
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 205755 68585 4.5758 0.01095 *
## Residuals 25 374713 14989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(ambient.MHV)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 205755 68585 4.5758 999 0.007 **
## Residuals 25 374713 14989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ambient3tr <- ambient %>%
mutate(logAlt = log(alt),
logbdo = log(bdo)) %>%
select(logAlt, oxy, logbdo)
ambient.pars3.d1 <- dist(ambient3tr)
(ambient.MHV2 <- betadisper(ambient.pars3.d1, gr))
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = ambient.pars3.d1, group = gr)
##
## No. of Positive Eigenvalues: 3
## No. of Negative Eigenvalues: 0
##
## Average distance to median:
## 1 2 3 4
## 8.016 8.984 9.630 7.336
##
## Eigenvalues for PCoA axes:
## PCoA1 PCoA2 PCoA3
## 13647.801 6.779 2.604
permutest(ambient.MHV2)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 17.44 5.814 0.0882 999 0.977
## Residuals 25 1647.78 65.911
Wilks.test(ambient3tr, gr)
##
## One-way MANOVA (Bartlett Chi2)
##
## data: x
## Wilks' Lambda = 0.077578, Chi2-Value = 62.633, DF = 9.000, p-value = 4.155e-10
## sample estimates:
## logAlt oxy logbdo
## 1 6.464881 113.88889 3.351171
## 2 6.245151 99.44444 3.512370
## 3 5.385688 83.87500 3.859659
## 4 5.477515 52.00000 5.010015
lw <- manova(as.matrix.data.frame(ambient3tr) ~ as.factor(gr))
summary(lw, test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## as.factor(gr) 3 0.077578 11.592 9 56.127 4.911e-10 ***
## Residuals 25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(peces.lda <- lda(gr ~ logAlt + oxy + logbdo, data = ambient3tr))
## Call:
## lda(gr ~ logAlt + oxy + logbdo, data = ambient3tr)
##
## Prior probabilities of groups:
## 1 2 3 4
## 0.3103448 0.3103448 0.2758621 0.1034483
##
## Group means:
## logAlt oxy logbdo
## 1 6.464881 113.88889 3.351171
## 2 6.245151 99.44444 3.512370
## 3 5.385688 83.87500 3.859659
## 4 5.477515 52.00000 5.010015
##
## Coefficients of linear discriminants:
## LD1 LD2 LD3
## logAlt -2.5129208 -1.72245625 -1.07719223
## oxy -0.0773069 -0.02117976 0.08197343
## logbdo -0.3348384 -2.71474681 2.30128498
##
## Proportion of trace:
## LD1 LD2 LD3
## 0.9032 0.0876 0.0092
peces.lda$means
## logAlt oxy logbdo
## 1 6.464881 113.88889 3.351171
## 2 6.245151 99.44444 3.512370
## 3 5.385688 83.87500 3.859659
## 4 5.477515 52.00000 5.010015
(C1 <- peces.lda$scaling)
## LD1 LD2 LD3
## logAlt -2.5129208 -1.72245625 -1.07719223
## oxy -0.0773069 -0.02117976 0.08197343
## logbdo -0.3348384 -2.71474681 2.30128498
nuevos_sitios <- data.frame(
logAlt = c(6.8, 5.5),
oxy = c(9, 10),
logbdo = c(3.2, 4.5)
)
nuevos_sitios
## logAlt oxy logbdo
## 1 6.8 9 3.2
## 2 5.5 10 4.5
(prediccion <- predict(peces.lda, newdata = nuevos_sitios))
## $class
## [1] 3 4
## Levels: 1 2 3 4
##
## $posterior
## 1 2 3 4
## 1 3.038541e-12 3.256050e-05 0.967908130 0.03205931
## 2 2.957066e-19 1.485777e-11 0.008073932 0.99192607
##
## $x
## LD1 LD2 LD3
## 1 4.780585 1.825078 -9.072993
## 2 7.534786 0.513920 -4.598999
ambient3trEstd <- as.data.frame(scale(ambient3tr))
peces.lda2 <- lda(gr ~ ., data = ambient3trEstd)
peces.lda2$means
## logAlt oxy logbdo
## 1 0.8194911 0.8683807 -0.5979420
## 2 0.4347680 0.2138834 -0.3316072
## 3 -1.0700611 -0.4915890 0.2421866
## 4 -0.9092809 -1.9358882 2.1428169
(C2 <- peces.lda2$scaling)
## LD1 LD2 LD3
## logAlt -1.4352215 -0.9837581 -0.6152241
## oxy -1.7061268 -0.4674274 1.8091148
## logbdo -0.2026608 -1.6430999 1.3928522
peces.lda2$svd^2
## [1] 53.6898562 5.2088449 0.5478896
(Fp2 <- predict(peces.lda2)$x)
## LD1 LD2 LD3
## 1 -4.08638573 -0.89640510 0.368028660
## 2 -2.49450633 0.46365901 -1.995824032
## 3 -2.80466835 -1.20357224 -0.404993616
## 4 -2.68895361 1.49616437 -2.201175466
## 5 -0.87807014 -2.09926524 -1.059020040
## 6 -2.51740981 -2.13333526 0.387269196
## 7 -2.90386963 0.07319673 -0.891988271
## 8 0.10415477 -1.24335518 -1.988893958
## 9 -1.49957974 -0.97965055 0.082158618
## 10 -1.88806718 0.38774388 0.504579631
## 11 -2.43308232 -0.02501061 1.334323268
## 12 -2.36655395 0.63877375 1.047520033
## 13 -2.35214071 -0.52520225 2.062059098
## 14 -1.57722535 1.28900168 0.253631503
## 15 -0.32438764 1.07783858 -0.206474231
## 16 -0.23770978 -0.21870301 1.018179031
## 17 -0.03051363 1.18888940 0.008410463
## 18 -0.14515674 0.79740513 0.706294065
## 19 0.34427132 1.44578197 0.169066314
## 20 1.44181530 0.83677120 0.075460197
## 21 1.38965445 0.44108254 0.553586732
## 22 3.22326372 -2.24627627 1.120313152
## 23 4.22156847 -1.19694492 -0.421313261
## 24 5.07604329 -1.72116640 -0.573615673
## 25 3.85542329 -0.32572787 -0.218162145
## 26 3.29378337 0.46604929 -0.152484763
## 27 2.84858565 1.28339090 -0.129929770
## 28 2.33552915 1.38947026 0.517474958
## 29 3.09418785 1.53939626 0.035520310
(peces.class2 <- predict(peces.lda2)$class)
## [1] 1 2 1 2 2 1 1 2 2 1 1 1 1 2 2 2 2 2 2 3 3 4 4 4 3 3 3 3 3
## Levels: 1 2 3 4
(peces.post2 <- round(predict(peces.lda2)$posterior, 2))
## 1 2 3 4
## 1 0.99 0.01 0.00 0.00
## 2 0.49 0.51 0.00 0.00
## 3 0.86 0.14 0.00 0.00
## 4 0.48 0.52 0.00 0.00
## 5 0.19 0.81 0.00 0.00
## 6 0.88 0.12 0.00 0.00
## 7 0.80 0.20 0.00 0.00
## 8 0.02 0.97 0.01 0.00
## 9 0.48 0.52 0.00 0.00
## 10 0.61 0.39 0.00 0.00
## 11 0.87 0.13 0.00 0.00
## 12 0.81 0.19 0.00 0.00
## 13 0.91 0.09 0.00 0.00
## 14 0.40 0.60 0.00 0.00
## 15 0.06 0.89 0.05 0.00
## 16 0.14 0.83 0.03 0.00
## 17 0.04 0.82 0.14 0.00
## 18 0.08 0.83 0.09 0.00
## 19 0.02 0.56 0.42 0.00
## 20 0.00 0.06 0.94 0.00
## 21 0.00 0.08 0.92 0.00
## 22 0.00 0.00 0.02 0.98
## 23 0.00 0.00 0.05 0.95
## 24 0.00 0.00 0.00 1.00
## 25 0.00 0.00 0.52 0.48
## 26 0.00 0.00 0.96 0.04
## 27 0.00 0.00 1.00 0.00
## 28 0.00 0.00 1.00 0.00
## 29 0.00 0.00 1.00 0.00
(peces.table2 <- table(gr, peces.class2))
## peces.class2
## gr 1 2 3 4
## 1 7 2 0 0
## 2 1 8 0 0
## 3 0 1 7 0
## 4 0 0 0 3
diag(prop.table(peces.table2, 1))
## 1 2 3 4
## 0.7777778 0.8888889 0.8750000 1.0000000
plot.lda(lda.out = peces.lda2,
groups = gr,
plot.sites = 2,
plot.centroids = 1,
mul.coef = 1.35
)
(peces.lda.jac <-
lda(gr ~ logAlt + oxy + logbdo,
data = ambient3trEstd,
CV = TRUE))
## $class
## [1] 1 2 1 1 2 1 1 2 2 1 1 1 1 1 2 2 2 2 2 3 3 4 4 4 4 3 3 3 3
## Levels: 1 2 3 4
##
## $posterior
## 1 2 3 4
## 1 9.856933e-01 1.430672e-02 1.672950e-10 2.340456e-16
## 2 2.052142e-01 7.947717e-01 1.408201e-05 1.575136e-11
## 3 8.264727e-01 1.735268e-01 4.779405e-07 3.687475e-11
## 4 9.296063e-01 7.039361e-02 1.335603e-07 1.981374e-18
## 5 3.608780e-01 6.387156e-01 3.810812e-04 2.530312e-05
## 6 9.921200e-01 7.879994e-03 1.463770e-09 1.358136e-11
## 7 7.506621e-01 2.493366e-01 1.326711e-06 2.232935e-12
## 8 2.540431e-02 9.231870e-01 5.048968e-02 9.190349e-04
## 9 4.093111e-01 5.905328e-01 1.560876e-04 2.306270e-08
## 10 5.758167e-01 4.240631e-01 1.202552e-04 6.768159e-10
## 11 8.452299e-01 1.547597e-01 1.041518e-05 5.549774e-11
## 12 7.725664e-01 2.274049e-01 2.867153e-05 3.153028e-11
## 13 8.673186e-01 1.326691e-01 1.224484e-05 2.925170e-10
## 14 5.108851e-01 4.880450e-01 1.069898e-03 3.173138e-10
## 15 7.462263e-02 8.627123e-01 6.266402e-02 1.043786e-06
## 16 1.894394e-01 7.565910e-01 5.394809e-02 2.154025e-05
## 17 4.593103e-02 7.585763e-01 1.954884e-01 4.215257e-06
## 18 9.989269e-02 7.611944e-01 1.389070e-01 5.986493e-06
## 19 1.464167e-02 8.008569e-01 1.845002e-01 1.212599e-06
## 20 3.740059e-04 7.139859e-02 9.276098e-01 6.175906e-04
## 21 7.829162e-04 1.031676e-01 8.944178e-01 1.631714e-03
## 22 3.852857e-07 3.730080e-04 1.441493e-01 8.554773e-01
## 23 7.456026e-10 7.136050e-06 8.107595e-02 9.189169e-01
## 24 4.629964e-13 3.321420e-08 6.661268e-03 9.933387e-01
## 25 1.217193e-09 2.224970e-05 2.948547e-01 7.051231e-01
## 26 7.089511e-08 2.921773e-04 9.501719e-01 4.953582e-02
## 27 3.279002e-07 6.627668e-04 9.965100e-01 2.826906e-03
## 28 4.432135e-06 2.499065e-03 9.967172e-01 7.793502e-04
## 29 5.892951e-08 2.011126e-04 9.973493e-01 2.449511e-03
##
## $terms
## gr ~ logAlt + oxy + logbdo
## attr(,"variables")
## list(gr, logAlt, oxy, logbdo)
## attr(,"factors")
## logAlt oxy logbdo
## gr 0 0 0
## logAlt 1 0 0
## oxy 0 1 0
## logbdo 0 0 1
## attr(,"term.labels")
## [1] "logAlt" "oxy" "logbdo"
## attr(,"order")
## [1] 1 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(gr, logAlt, oxy, logbdo)
## attr(,"dataClasses")
## gr logAlt oxy logbdo
## "numeric" "numeric" "numeric" "numeric"
##
## $call
## lda(formula = gr ~ logAlt + oxy + logbdo, data = ambient3trEstd,
## CV = TRUE)
##
## $xlevels
## named list()
peces.jac.class <- peces.lda.jac$class
peces.jac.table <- table(gr, peces.jac.class)
diag(prop.table(peces.jac.table, 1))
## 1 2 3 4
## 0.7777778 0.6666667 0.7500000 1.0000000