Practica Grad Indirecto
Federico J. Villatoro
2023-09-14
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")}
# 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)
require(stats)
install_github("vqv/ggbiplot")
require(ggbiplot)Data
data(dune)
data(dune.env)PCA
M1 <- dune
dune.pca <- prcomp(M1, center = TRUE, scale. = TRUE)
fviz_pca_ind(dune.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = dune.env$Management, # color by groups
palette = c("#00AFBB", "#E7B800", "#FC4E07", "#E7B880"),
addEllipses = FALSE, # Concentration ellipses
legend.title = "Management"
)
fviz_pca_ind(dune.pca, geom.ind = "point", col.ind = dune.env$Management,
palette = c("#00AFBB", "#E7B800", "#FC4E07", "#E7B880"),
addEllipses = TRUE, ellipse.type = "confidence",
legend.title = "Management"
)
fviz_pca_var(dune.pca, col.var = "red", repel = TRUE)
ind.p <- fviz_pca_ind(dune.pca, geom = "point", col.ind = dune.env$Management)
ggpubr::ggpar(ind.p,
title = "Principal Component Analysis",
subtitle = "Dune data set",
caption = "Source: factoextra",
xlab = "PC1", ylab = "PC2",
legend.title = "Management", legend.position = "top",
ggtheme = theme_gray(), palette = "jco"
)
fviz_pca_biplot(dune.pca,
col.ind = dune.env$Management, palette = "jco",
addEllipses = TRUE,
ellipse.type= "confidence",
label = "var",
col.var = "black", repel = TRUE,
legend.title = "Management") 
fviz_pca_biplot(dune.pca,
# Fill individuals by groups
geom.ind = "point",
pointshape = 21,
pointsize = 2.5,
fill.ind = dune.env$Management,
col.ind = "black",
# Color variable by groups
#col.var = ,
legend.title = list(fill = "Management", color = "Clusters"),
repel = TRUE # Avoid label overplotting
)+
ggpubr::fill_palette("jco")+ # Indiviual fill color
ggpubr::color_palette("npg") # Variable colors
### Correlations with PCs
DUNE <- dune
DUNE$PCA1 <- dune.pca$x[,1]
DUNE$PCA2 <- dune.pca$x[,2]
DUNE$PCA3 <- dune.pca$x[,3]
DUNE$PCA4 <- dune.pca$x[,4]
print(DUNE[,c(31:34)])## PCA1 PCA2 PCA3 PCA4
## 1 0.5276579 1.08011979 0.0702171 -0.77828260
## 2 3.1866021 1.74654582 0.4008214 -2.29580010
## 3 0.6437881 2.47303132 1.0262201 -0.71220528
## 4 0.4974602 3.17153838 2.9214282 -2.12522365
## 5 3.3440281 -0.04613685 -1.5171195 1.21500338
## 6 3.4505427 -1.08270271 -3.1738294 2.99981695
## 7 2.8511292 0.07723272 -1.4205363 1.48317091
## 8 -1.8502140 1.21148541 0.1548990 -0.14693872
## 9 -0.3839450 2.29309029 1.1941347 1.50056566
## 10 3.9497450 -0.16503613 -1.1943440 -1.90361222
## 11 1.5579455 -1.67360338 0.3259882 -1.45114800
## 12 -1.2481921 1.32881704 1.7519094 2.67363043
## 13 -1.6317972 2.52839510 2.0415588 2.61111669
## 14 -3.1861923 -0.51959100 -2.1601829 -1.11188282
## 15 -3.4547894 -0.64823033 -1.7277829 -0.52734223
## 16 -4.5543039 0.54810093 -1.6006539 -0.03556033
## 17 0.3823832 -2.71434036 0.6823393 0.41084792
## 18 0.7903125 -1.75439962 -0.4725707 -1.45372867
## 19 -0.1982217 -6.51396167 4.4010946 0.46907513
## 20 -4.6739390 -1.34035476 -1.7035915 -0.82150246#
PCA1.cor <- cor(DUNE[,c(1:30)],DUNE$PCA1)
PCA2.cor <- cor(DUNE[,c(1:30)],DUNE$PCA2)
PCA3.cor <- cor(DUNE[,c(1:30)],DUNE$PCA3)
PCA4.cor <- cor(DUNE[,c(1:30)],DUNE$PCA4)
sp.cor <- as.data.frame(cbind(PCA1.cor,
PCA2.cor,
PCA3.cor,
PCA4.cor))
sp.cor## V1 V2 V3 V4
## Achimill 0.735519493 -0.035557766 -0.323423120 -0.079958862
## Agrostol -0.717190544 0.503736312 0.141687250 -0.037393221
## Airaprae 0.004292982 -0.747591155 0.517127489 0.091741180
## Alopgeni -0.263514202 0.587291705 0.382488500 0.298599002
## Anthodor 0.536437669 -0.565560825 -0.047187096 0.242308636
## Bellpere 0.534545267 0.285192236 0.073269896 -0.524328510
## Bromhord 0.597399534 0.265626299 -0.005646541 -0.408213944
## Chenalbu -0.144835270 0.266218262 0.254869483 0.377966241
## Cirsarve 0.044153634 0.333935718 0.364712932 -0.307631902
## Comapalu -0.428219040 -0.089329491 -0.352617549 -0.172381582
## Eleopalu -0.778583238 -0.014472318 -0.429279957 -0.154244912
## Elymrepe 0.290974684 0.537200973 0.251924873 -0.152232947
## Empenigr -0.017593785 -0.685864148 0.549435409 0.067899900
## Hyporadi 0.046459631 -0.788607181 0.543256509 0.006942918
## Juncarti -0.633296249 0.121524606 -0.196626850 0.008792643
## Juncbufo -0.081874742 0.377309150 0.303265733 0.641684995
## Lolipere 0.691559075 0.283887914 -0.100357788 -0.281382071
## Planlanc 0.691177714 -0.259258999 -0.475653050 0.247174326
## Poaprat 0.656568603 0.472589572 0.086667046 -0.240620140
## Poatriv 0.391483428 0.716617533 0.186469900 0.317394766
## Ranuflam -0.812111853 0.017621551 -0.318514659 -0.047057469
## Rumeacet 0.470035653 0.009772765 -0.386928223 0.661742324
## Sagiproc -0.098421938 0.156605916 0.758015566 0.135273074
## Salirepe -0.307315312 -0.531927662 0.065409574 -0.163899116
## Scorautu 0.453235138 -0.518811796 0.255584724 -0.175410374
## Trifprat 0.476615569 -0.102068443 -0.491747383 0.534807128
## Trifrepe 0.443919007 0.043766766 -0.213912989 -0.036216733
## Vicilath 0.298034750 -0.237021027 -0.054157137 -0.387289750
## Bracruta -0.068122718 -0.317322531 -0.266146556 0.147101858
## Callcusp -0.649485246 -0.084901555 -0.419221776 -0.183885477sp.cor[order(sp.cor$V1),] #Orden de acuerdo a los valores de los coeficientes de correlación con el primer componente (V1 en este data.frame)## V1 V2 V3 V4
## Ranuflam -0.812111853 0.017621551 -0.318514659 -0.047057469
## Eleopalu -0.778583238 -0.014472318 -0.429279957 -0.154244912
## Agrostol -0.717190544 0.503736312 0.141687250 -0.037393221
## Callcusp -0.649485246 -0.084901555 -0.419221776 -0.183885477
## Juncarti -0.633296249 0.121524606 -0.196626850 0.008792643
## Comapalu -0.428219040 -0.089329491 -0.352617549 -0.172381582
## Salirepe -0.307315312 -0.531927662 0.065409574 -0.163899116
## Alopgeni -0.263514202 0.587291705 0.382488500 0.298599002
## Chenalbu -0.144835270 0.266218262 0.254869483 0.377966241
## Sagiproc -0.098421938 0.156605916 0.758015566 0.135273074
## Juncbufo -0.081874742 0.377309150 0.303265733 0.641684995
## Bracruta -0.068122718 -0.317322531 -0.266146556 0.147101858
## Empenigr -0.017593785 -0.685864148 0.549435409 0.067899900
## Airaprae 0.004292982 -0.747591155 0.517127489 0.091741180
## Cirsarve 0.044153634 0.333935718 0.364712932 -0.307631902
## Hyporadi 0.046459631 -0.788607181 0.543256509 0.006942918
## Elymrepe 0.290974684 0.537200973 0.251924873 -0.152232947
## Vicilath 0.298034750 -0.237021027 -0.054157137 -0.387289750
## Poatriv 0.391483428 0.716617533 0.186469900 0.317394766
## Trifrepe 0.443919007 0.043766766 -0.213912989 -0.036216733
## Scorautu 0.453235138 -0.518811796 0.255584724 -0.175410374
## Rumeacet 0.470035653 0.009772765 -0.386928223 0.661742324
## Trifprat 0.476615569 -0.102068443 -0.491747383 0.534807128
## Bellpere 0.534545267 0.285192236 0.073269896 -0.524328510
## Anthodor 0.536437669 -0.565560825 -0.047187096 0.242308636
## Bromhord 0.597399534 0.265626299 -0.005646541 -0.408213944
## Poaprat 0.656568603 0.472589572 0.086667046 -0.240620140
## Planlanc 0.691177714 -0.259258999 -0.475653050 0.247174326
## Lolipere 0.691559075 0.283887914 -0.100357788 -0.281382071
## Achimill 0.735519493 -0.035557766 -0.323423120 -0.079958862Graficar Cuatro especies por Componente según los datos de correlación calculados arriba. Con el primer eje (PCA1), las primeras dos especies ya están incluidas. Continuar escribiendo los nombres de las especies que se correlacionan más (cor > 0.7), con el PCA1, y luego con el PCA2, y el PCA3 y el PCA4.
# PCA1
p1 <- ggplot(DUNE, aes(x= PCA1, y = Ranuflam ))
p1 <- p1 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p2 <- ggplot(DUNE, aes(x= PCA1, y = Achimill ))
p2 <- p2 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p3 <- ggplot(DUNE, aes(x= PCA1, y = ))
p3 <- p3 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p4 <- ggplot(DUNE, aes(x= PCA1, y = ))
p4 <- p4 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
grid.arrange(p1,p2,p3,p4)
# PCA2
p1 <- ggplot(DUNE, aes(x= PCA2, y = ))
p1 <- p1 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p2 <- ggplot(DUNE, aes(x= PCA2, y = ))
p2 <- p2 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p3 <- ggplot(DUNE, aes(x= PCA2, y = ))
p3 <- p3 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p4 <- ggplot(DUNE, aes(x= PCA2, y = ))
p4 <- p4 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
grid.arrange(p1,p2,p3,p4)
# PCA3
p1 <- ggplot(DUNE, aes(x= PCA3, y = ))
p1 <- p1 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p2 <- ggplot(DUNE, aes(x= PCA3, y = ))
p2 <- p2 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p3 <- ggplot(DUNE, aes(x= PCA3, y = ))
p3 <- p3 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p4 <- ggplot(DUNE, aes(x= PCA3, y = ))
p4 <- p4 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
grid.arrange(p1,p2,p3,p4)
# PCA4
p1 <- ggplot(DUNE, aes(x= PCA4, y = ))
p1 <- p1 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p2 <- ggplot(DUNE, aes(x= PCA4, y = ))
p2 <- p2 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p3 <- ggplot(DUNE, aes(x= PCA4, y = ))
p3 <- p3 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
p4 <- ggplot(DUNE, aes(x= PCA4, y = ))
p4 <- p4 + geom_point(position = position_jitter(width= 0.2, height= 0.2 ), size=1 ) +
geom_smooth(method="glm", method.args=list(family="poisson"), se=FALSE)
grid.arrange(p1,p2,p3,p4)