library(FactoMineR)
library(ggplot2)
library(factoextra)

data = data # built in data set
head(data) # 4 continuous variables and 1 categorical variable



data2 <- data[1:18]
head(data2)


# PCA analysis to get PCs
data.pca <- PCA(data2, scale.unit = TRUE, graph = FALSE)


# look at eigenvalues (measure of how much variance is contained in a PC)
# there are the same number of PCs as original variables
data.pca$eig


# Scree plot with factoextra
fviz_screeplot(data.pca, ncp = 18) # ncp = number of PCs to show

# Simple PCA factor map with FactoMineR graphics
plot.PCA(data.pca, axes = c(1,2), choix = "var")


# Biplot with factoextra
# Biplot means it has individuals and contributions
fviz_pca(data.pca)


# include only contributions of variables to clean it up
fviz_pca_var(data.pca, col.var = "contrib")


# Control scale colors
fviz_pca_var(data.pca, col.var = "contrib") +
  scale_color_gradient2(low = "blue", mid = "steelblue", high = "red", midpoint = 25.0) +
  theme_void()

# include only individuals with no labels
fviz_pca_ind(data.pca, label = "none")

# use Species from data to change habillage
fviz_pca_ind(data.pca, label="none", habillage = as.factor(data$Sites))
#install.packages("ggforce")
library(ggforce)
# changed color and shape

# add ellipses
fviz_pca_ind(data.pca, label="none", habillage = as.factor(data$Sites),
             repel = TRUE, 
             # Don't use default Ellipses!!!!
             # addEllipses = TRUE,
             invisible='quali') +
  # ADD ggforce's ellipses
  ggforce::geom_mark_ellipse(aes(fill = Groups,
                                 color = Groups)) +
  theme(legend.position = 'right') +
  coord_equal()



# Make a pretty biplot
fviz_pca_biplot(data.pca,
                # individuals
                geom.ind = "point",
                fill.ind = data$Sites, col.ind = "black",
                pointshape = 21, pointsize = 2,
                addEllipses = TRUE, repel= TRUE,
                # variables
                col.var = "contrib",
                gradient.cols = "RdYlBu",
                legend.title = list(fill = "Sites",
                                    color = "Contrib",
                                    alpha = "Contrib")) +
  scale_color_gradient2(low = "red", mid = "steelblue", high = "red", midpoint = 25.0)



## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.


A <- fviz_pca_biplot(data.pca,
                     # individuals
                     geom.ind = "point",
                     fill.ind = data$Sites, col.ind = "black",
                     pointshape = 21, pointsize = 2,
                     addEllipses = TRUE, repel= TRUE,
                     # variables
                     col.var = "contrib",
                     gradient.cols = "RdYlBu",
                     legend.title = list(fill = "Sites",
                                         color = "Contrib",
                                         alpha = "Contrib")) +
  scale_color_gradient2(low = "red", mid = "steelblue", high = "red", midpoint = 25.0)+
  theme_minimal()+
   theme(legend.position = "right")+
  ggtitle("")+
  xlab("PC1 (16.7%)")+
  ylab("PC2 (14.5%)")+
  theme(axis.title.x=element_text(size=8, color="black", face="bold"),
        axis.title.y=element_text(size=8, color="black", face="bold"),
        plot.title=element_text(size=10, face="bold", color="red"))




ggsave("Water and Air Quality.pdf", A, limitsize = FALSE)