Trabajo_Final
Federico J. Villatoro
2023-11-07
Instrucciones: Correr el código contenido en este Rpubs, completando lo que haga falta para obtener el mejor modelo de variables ambientales que explique la matriz de especies.
packages
# install if required
if(!require(tidyverse)){install.packages("tidyverse")}
if(!require(googlesheets4)){install.packages("googlesheets4")}
if(!require(GGally)){install.packages("GGally")}
if(!require(reshape2)){install.packages("reshape2")}
if(!require(lme4)){install.packages("lme4")}
if(!require(compiler)){install.packages("compiler")}
if(!require(parallel)){install.packages("parallel")}
if(!require(boot)){install.packages("boot")}
if(!require(lattice)){install.packages("lattice")}
if(!require(car)){install.packages("car")}
if(!require(ggpubr)){install.packages("ggpubr")}
if(!require(vegan)){install.packages("vegan")}
if(!require(ggsci)){install.packages("ggsci")}
if(!require(ggrepel)){install.packages("ggrepel")}
if(!require(ggforce)){install.packages("ggforce")}
if(!require(plotly)){install.packages("plotly")}
if(!require(rgl)){install.packages("rgl")}
if(!require(mgcv)){install.packages("mgcv")}
if(!require(factoextra)){install.packages("factoextra")}
if(!require(devtools)){install.packages("devtools")}
if(!require(gridExtra)){install.packages("gridExtra")}
if(!require(BiodiversityR)){install.packages("BiodiversityR")}
# Load
library(tidyverse)
library(googlesheets4); gs4_deauth()
library(googledrive)
library(GGally)
library(reshape2)
library(lme4)
library(compiler)
library(parallel)
library(boot)
library(lattice)
library(car)
library(ggpubr)
library(vegan)
library(readxl)
library(ggsci)
library(ggrepel)
library(ggforce)
library(plotly)
library(rgl)
library(mgcv)
library(dplyr)
library(factoextra)
library(devtools)
library(gridExtra)
library(BiodiversityR) # also loads vegan
library(readxl)
library(ggsci)
library(ggrepel)
library(ggforce)
require(stats)
install_github("vqv/ggbiplot")
require(ggbiplot)Dataset
Datos asociados al estudio de Condit et al. 2002. Beta-diversity in tropical forest trees. Science 295: 666-669. Supplementary material [data] (https://www.davidzeleny.net/anadat-r/doku.php/en:data:bci)
data(BCI)
BCI.env <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/bci.env.txt', row.names = 1)
BCI.soil <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/bci.soil.txt')Matrices
Se quiere evaluar la combinación de variables (ambientales u de química de suelo) que mejor explique la diversidad de especies vegetales en este paisaje.
M1 <- BCI
M2 <- cbind(BCI.env[,3:5], BCI.soil[,c(4,8,15)])
# checking, rownames
str(M2)## 'data.frame': 50 obs. of 6 variables:
## $ elevation: num 130 137 144 147 144 ...
## $ convex : num -7.87 -10.7 -14.67 -16.76 -12.48 ...
## $ slope : num 6.69 5.09 3.1 1.87 5.12 ...
## $ B : num 0.794 0.67 0.595 0.568 0.399 ...
## $ K : num 141.9 137.2 98.7 98.4 94.1 ...
## $ pH : num 4.32 4.38 4.35 4.46 4.4 ...check.datasets(M1, M2)## OKvegan
UN-CONSTRAINED
Chequeamos el largo del gradiente ambiental al que las especies están respondiendo. Gradientes > 3 se consideran “grandes” y se prefiere utilizar un método de ordenación que asuma respuestas unimodales de las especies. Es decir, que utilice promedios ponderados para calcular los puntajes de los sitios (e.g., ordenación unimodal), y no sumas ponderadas (e.g., ordenación lineal). Utilizar aquí el método que nos sirve para calcular el largo de los gradientes ambientales (y la biodiversidad beta potencial que logramos captar).
Basándonos en el largo del gradiente, Utilizar ya sea:
- CCA
ó - RDA
Ordination Model
MODEL Building [Vegan PCKG functions]
upr <- vegan::____(M1~ ., data=M2) # Full Model
lwr <- vegan::____(M1~ 1, data=M2) # Null model
vif.cca(upr) # V.I.F for the full modelEste es el método más recomendado según Blanchet, Legendre & Borcard. 2008. Ecology 89, 2623-2623 para seleccionar las variables del modelo.
Es un ‘Forward selection’ pero se detiene una vez el r^2 del modelo globas es sobrepasado por el modelo que se está poniendo a prueba.
## Stepwise Selection via Adjusted R^2
### 1. Global test
### 2. Radjusted forward selection
finalModel <- ordiR2step(lwr, scope= formula(upr), trace = FALSE)
finalModel$anova
# Adjusted R^2
RsquareAdj(upr)
RsquareAdj(finalModel)
# VIF
vif.cca(finalModel) # verificar si estamos "OK" con el Factor de inflación de la varianza ?Permutest
Este método para probar significancia de cada témrino, es el menos sesgado.
Calcula la importancia de cada término, cuando este se remueve del modelo global, i.e., “Efecto Marginal”
### Marginal effects (MORE USEFUL !)
set.seed(10)
anova(finalModel, by = "margin")Diagrama de ordenación
.long formats
plot2 <- ordiplot(finalModel, choices = c(1,2), scaling="none")
sites.long <- sites.long(plot2, env.data=M2)
species.long <- species.long(plot2)
axis.long <- axis.long(finalModel, choices = c(1,2))env.fit (Spp)
spec.envfit <- envfit(plot2, env=M1)
spec.data.envfit <- data.frame(r=spec.envfit$vectors$r, p=spec.envfit$vectors$pvals)
species.long2 <- species.long(plot2, spec.data=spec.data.envfit)
species.long2En este caso no utilizamos envfit() para las
variables ambientales. Necesitamos, en cambio, la obtención de los
‘biplot scores’ resultantes del análisis de ordenación. En el
código siguiente se obtienen dichos ‘biplot scores’.
## Vector Fit
vectors<- data.frame(cbind(vector= rownames(data.frame(plot2$biplot)), axis1=data.frame(plot2$biplot)[,1] , axis2=data.frame(plot2$biplot)[,2] ) )
vectors$vector <- sapply(vectors[,1], as.factor)
vectors$axis1 <- as.numeric(vectors$axis1)
vectors$axis2 <- as.numeric(vectors$axis2)
vectors
#vectors.envfit <- envfit(plot2, env= M2[,c(3,14)])
#vectors.long3 <- vectorfit.long()
#vectors.long3### Spec r >= 30%
species.long3 <- species.long2[species.long2$r >= 0.3, ]
species.long3
### Var Ambientales r >= 0
#vectors.long3 <- vectors.long3[vectors.long3$r >= 0, ]
#vectors.long3ggplots
f = 1
plotgg2 <- ggplot() +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
xlab(axis.long[1, "label"]) +
ylab(axis.long[2, "label"]) +
scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_point(data=sites.long,
aes(x=axis1, y=axis2),
size=3, shape=1) +
#geom_text_repel(data=sites.long,aes(x=axis1, y=axis2, label=labels),colour="black") +
## Spp
geom_point(data=species.long3,aes(x=axis1, y=axis2)) +
geom_segment(data=species.long3,aes(x=0, y=0, xend=axis1, yend=axis2),colour="red", linewidth=0.7, arrow=arrow(angle=20)) +
geom_text_repel(data=species.long3,aes(x=axis1, y=axis2, label=labels),colour="red") +
## ENV
geom_segment(data=vectors,aes(x=0, y=0, xend=axis1*f, yend=axis2*f),colour="blue", linewidth=0.7, arrow=arrow(angle=20)) +
geom_text_repel(data=vectors,aes(x=axis1*f, y= axis2*f, label=vector),colour="blue") +
ggsci::scale_colour_npg() +
coord_fixed(ratio=1)
plotgg2Estas figura que se muestra aquí está ocultando los nombres de las variables que resultan del modelo. Pero el código , al ejecutarlo, sí mostrará los nombres en el plot (también para los ordisurfs).
Ordisurfs
plot2 <- ordiplot(finalModel, choices = c(1,2), display = "sites", scaling=0)
ordisurf(plot2, y=as.numeric(M2$elevation), add=TRUE, col="blue", main= ) # "Elevation" species response
## species response
ordisurf(plot2, y=as.numeric(M1$Trichilia.tuberculata), add=TRUE, col="orange", main= "'Trichilia.tuberculata (orange)' response to elevation gradient") # "Trichilia.tuberculata" species response
plot2 <- ordiplot(finalModel, choices = c(1,2), display = "sites", scaling=0)
ordisurf(plot2, y=as.numeric(M2$slope), add=TRUE, col="blue", main= ) # "Elevation" species response
## species response
ordisurf(plot2, y=as.numeric(M1$Gustavia.superba), add=TRUE, col="red", main= "'Gustavia.superba (red)' response to slope gradient") # "Gustavia.superba" species response