library(dplR)
library(ggplot2)
library(pointRes)
library(reshape2)
library(detrendeR) 
library(pointRes) 

setwd("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA")

###1.Chronologies ----

crono.TNFLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFLN.rwl")
crono.TNFIP<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFIP.rwl")
crono.TNFPE<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFPE.rwl")
crono.TNFPL<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFPL.rwl")
crono.TNFPI<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFPI.rwl")
crono.LAPLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAP/LAPLNd.rwl")
crono.LAPPI<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAP/LAPPId.rwl")
crono.LAGLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAG/LAGLNd.rwl")
crono.LAGPI<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAG/LAGPId.rwl")
crono.AZLA<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/AZO/AZLA.rwl")
crono.AZIA<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/AZO/AZIA.rwl")
crono.AZJB<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/AZO/AZJB.rwl")
crono.MADLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/MAD/MADLNd.rwl")
crono.MADOF<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/MAD/MADOFd.rwl")
crono.MADCA<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/MAD/MADCAd.rwl")

###2.Residuals (RWI)----

#########___ISLAS CANARIAS----
#___TENERIFE----

crono.TNFIP<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFIP.rwl")
TNFIP.cro <- chron(crono.TNFIP)
plot(TNFIP.cro)
negexpcrono.TNFIP<-detrend(crono.TNFIP,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.TNFIP<-detrend(negexpcrono.TNFIP,method='Spline')
spline.crn.TNFIP<-chron(splinegxpcrono.TNFIP)
plot(spline.crn.TNFIP)

crono.TNFPE<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFPE.rwl")
TNFPE.cro <- chron(crono.TNFPE)
plot(TNFPE.cro)
negexpcrono.TNFPE<-detrend(crono.TNFPE,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.TNFPE<-detrend(negexpcrono.TNFPE,method='Spline')
spline.crn.TNFPE<-chron(splinegxpcrono.TNFPE)
plot(spline.crn.TNFPE) 

crono.TNFPL<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFPL.rwl")
TNFPL.cro <- chron(crono.TNFPL)
plot(TNFPL.cro)
negexpcrono.TNFPL<-detrend(crono.TNFPL,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.TNFPL<-detrend(negexpcrono.TNFPL,method='Spline')
spline.crn.TNFPL<-chron(splinegxpcrono.TNFPL)
plot(spline.crn.TNFPL)

crono.TNFPI<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFPI.rwl")
TNFPI.cro <- chron(crono.TNFPI)
plot(TNFPI.cro)
negexpcrono.TNFPI<-detrend(crono.TNFPI,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.TNFPI<-detrend(negexpcrono.TNFPI,method='Spline')
spline.crn.TNFPI<-chron(splinegxpcrono.TNFPI)
plot(spline.crn.TNFPI)

crono.TNFLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/TNF/TNFLN.rwl")
TNFLN.cro <- chron(crono.TNFLN)
plot(TNFLN.cro)
negexpcrono.TNFLN<-detrend(crono.TNFLN,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.TNFLN<-detrend(negexpcrono.TNFLN,method='Spline')
spline.crn.TNFLN<-chron(splinegxpcrono.TNFLN)
plot(spline.crn.TNFLN)

#___LA PALMA----

crono.LAPLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAP/LAPLNd.rwl")
LAPLN.cro <- chron(crono.LAPLN)
plot(LAPLN.cro)
negexpcrono.LAPLN<-detrend(crono.LAPLN,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.LAPLN<-detrend(negexpcrono.LAPLN,method='Spline')
spline.crn.LAPLN<-chron(splinegxpcrono.LAPLN)
plot(spline.crn.LAPLN)

crono.LAPPI<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAP/LAPPId.rwl")
LAPPI.cro <- chron(crono.LAPPI)
plot(LAPPI.cro)
negexpcrono.LAPPI<-detrend(crono.LAPPI,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.LAPPI<-detrend(negexpcrono.LAPPI,method='Spline')
spline.crn.LAPPI<-chron(splinegxpcrono.LAPPI)
plot(spline.crn.LAPPI)

#___LA GOMERA----

crono.LAGLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAG/LAGLNd.rwl")
LAGLN.cro <- chron(crono.LAGLN)
plot(LAGLN.cro)
negexpcrono.LAGLN<-detrend(crono.LAGLN,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.LAGLN<-detrend(negexpcrono.LAGLN,method='Spline')
spline.crn.LAGLN<-chron(splinegxpcrono.LAGLN)
plot(spline.crn.LAGLN)

crono.LAGPI<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/LAG/LAGPId.rwl")
LAGPI.cro <- chron(crono.LAGPI)
plot(LAGPI.cro)
negexpcrono.LAGPI<-detrend(crono.LAGPI,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.LAGPI<-detrend(negexpcrono.LAGPI,method='Spline')
spline.crn.LAGPI<-chron(splinegxpcrono.LAGPI)
plot(spline.crn.LAGPI)

#write table__CANARIAS----  

# write.table(spline.crn.TNFIP,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\TNFIP.STD.csv")
# write.table(spline.crn.TNFPE,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\TNFPE.STD.csv")
# write.table(spline.crn.TNFPL,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\TNFPL.STD.csv")
# write.table(spline.crn.TNFPI,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\TNFPI.STD.csv")
# write.table(spline.crn.TNFLN,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\TNFLN.STD.csv")
# write.table(spline.crn.LAPLN,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\LAPLN.STD.csv")
# write.table(spline.crn.LAPPI,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\LAPPI.STD.csv")
# write.table(spline.crn.LAGLN,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\LAGLN.STD.csv")
# write.table(spline.crn.LAGPI,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\LAGPI.STD.csv")


#########___ISLAS AZORES(Terceira)----

crono.AZLA<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/AZO/AZLA.rwl")
AZLA.cro <- chron(crono.AZLA)
plot(AZLA.cro)
negexpcrono.AZLA<-detrend(crono.AZLA,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.AZLA<-detrend(negexpcrono.AZLA,method='Spline')
spline.crn.AZLA<-chron(splinegxpcrono.AZLA)
plot(spline.crn.AZLA)

crono.AZIA<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/AZO/AZIA.rwl")
AZIA.cro <- chron(crono.AZIA)
plot(AZIA.cro)
negexpcrono.AZIA<-detrend(crono.AZIA,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.AZIA<-detrend(negexpcrono.AZIA,method='Spline')
spline.crn.AZIA<-chron(splinegxpcrono.AZIA)
plot(spline.crn.AZIA) 

crono.AZJB<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/AZO/AZJB.rwl")
AZJB.cro <- chron(crono.AZJB)
plot(AZJB.cro)

negexpcrono.AZJB<-detrend(crono.AZJB,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.AZJB<-detrend(negexpcrono.AZJB,method='Spline')
spline.crn.AZJB<-chron(splinegxpcrono.AZJB)
plot(spline.crn.AZJB)

#write table__AZORES----  

#write.table(spline.crn.AZLA,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\AZLA.STD.csv")
# write.table(spline.crn.AZIA,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\AZIA.STD.csv")
# write.table(spline.crn.AZJB,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\AZJB.STD.csv")

#########___MADEIRA----  

crono.MADLN<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/MAD/MADLNd.rwl")
MADLN.cro <- chron(crono.MADLN)
plot(MADLN.cro)
negexpcrono.MADLN<-detrend(crono.MADLN,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.MADLN<-detrend(negexpcrono.MADLN,method='Spline')
spline.crn.MADLN<-chron(splinegxpcrono.MADLN)
plot(spline.crn.MADLN)

crono.MADOF<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/MAD/MADOFd.rwl")
MADOF.cro <- chron(crono.MADOF)
plot(MADOF.cro)
negexpcrono.MADOF<-detrend(crono.MADOF,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.MADOF<-detrend(negexpcrono.MADOF,method='Spline')
spline.crn.MADOF<-chron(splinegxpcrono.MADOF)
plot(spline.crn.MADOF)    

crono.MADCA<-read.rwl("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/SERIES LAURISILVA RESPUESTA CLIMA/CRONOLOGIAS FINALES LAURISILVA/MAD/MADCAd.rwl")
MADCA.cro <- chron(crono.MADCA)
plot(MADCA.cro)
negexpcrono.MADCA<-detrend(crono.MADCA,method='ModNegExp',constrain.nls='when.fail')
splinegxpcrono.MADCA<-detrend(negexpcrono.MADCA,method='Spline')
spline.crn.MADCA<-chron(splinegxpcrono.MADCA)
plot(spline.crn.MADCA) 

#write table__MADEIRA----  

# write.table(spline.crn.MADLN,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\MADLN.STD.csv")
# write.table(spline.crn.MADOF,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\MADOF.STD.csv")
# write.table(spline.crn.MADCA,"C:\\Users\\Usuario\\Desktop\\ART?CULO RESPUESTA CLIMA\\SERIES LAURISILVA RESPUESTA CLIMA\\cronos.std\\MADCA.STD.csv")

###3. Statistics with detrendeR----

install.packages("detrendeR")
library(detrendeR)

detrender()

###4.Principal Component Analysis (PCA)----

library('corrr')
library(ggcorrplot)
library("FactoMineR")
library(factoextra)
library(psych)

#PCA_clima___AZORES----

pca.azo <- read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/graficos corr-PCA/PCA_LIMPIOS CON ORLSS TABLA/PCA.ter.csv")
corr.azo.cre <- as.data.frame(lapply(pca.azo, as.numeric))
str(corr.azo.cre)
corr_matrix <- cor(corr.azo.cre)
head(corr_matrix)
ggcorrplot(corr_matrix)
p_value_matrix<- corr.test(corr.azo.cre)$p #clacular matriz con p-valores
print(p_value_matrix)# Mostrar la matriz de p-valores
data.pca <- princomp(corr_matrix)
summary(data.pca)

data.pca$loadings[, 1:2] #cargar/mostrar componentes principales 1 y 2
fviz_eig(data.pca, addlabels = TRUE)

grp <- c("JB","LA","IA", "P.summer", "NAO.annual", "EA.winter")

fviz_pca_var(data.pca, col.var = grp, palette = c("darkblue","darkgreen", "darkgreen","darkgreen","darkblue","darkred"), labelsize = 5, repel = TRUE, arrowsize = 1) + 
  theme(panel.grid = element_blank(),text = element_text(size = 18), axis.title = element_text(size = 20),
        axis.text = element_text(size = 20)) 


#PCA_clima___MADEIRA----

pca.mad<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/graficos corr-PCA/PCA_LIMPIOS CON ORLSS TABLA/PCA.mad.csv")
corr.mad.cre <- as.data.frame(lapply(pca.mad, as.numeric))
str(corr.mad.cre)
corr_matrix <- cor(corr.mad.cre)
head(corr_matrix)
ggcorrplot(corr_matrix)
p_value_matrix<- corr.test(corr.mad.cre)$p #clacular matriz con p-valores
print(p_value_matrix)# Mostrar la matriz de p-valores
data.pca <- princomp(corr_matrix)
summary(data.pca)
data.pca$loadings[, 1:2] #cargar/mostrar componentes principales 1 y 2
fviz_eig(data.pca, addlabels = TRUE)


grp <- c("LN","CA","OF", "SPEI.winter", "P.winter", "P.PH", "NAO.HY", "Tme.autumn", "Tme.winter")

fviz_pca_var(data.pca, col.var = grp, palette = c("darkgreen","darkgreen", "darkred" ,"darkgreen", "darkred", "darkred","darkred","darkred", "darkblue"), labelsize = 5, repel = TRUE, arrowsize = 1) + 
  theme(panel.grid = element_blank(),text = element_text(size = 18), axis.title = element_text(size = 20),
        axis.text = element_text(size = 20))                                                                                 

#PCA_clima___TENERIFE----

pca.tnf <- read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/graficos corr-PCA/PCA_LIMPIOS CON ORLSS TABLA/PCA.tnf.csv")
corr.TNF.cre <- as.data.frame(lapply(pca.tnf, as.numeric))
str(corr.TNF.cre)
corr_matrix <- cor(corr.TNF.cre)
head(corr_matrix)
ggcorrplot(corr_matrix)
p_value_matrix<- corr.test(corr.TNF.cre)$p #clacular matriz con p-valores
print(p_value_matrix)# Mostrar la matriz de p-valores
data.pca <- princomp(corr_matrix)
summary(data.pca)
data.pca$loadings[, 1:2] #cargar/mostrar componentes principales 1 y 2
fviz_eig(data.pca, addlabels = TRUE)

grp <- c("LN","PI","IP", "PE", "PL", "SPEI.HY", "SPEI.winter", "P.winter", "EA.HY", "Tme.summer", "EA.autumn", "NAO.winter","NAO.spring")

fviz_pca_var(data.pca, col.var = grp, palette = c("darkblue","darkblue", "darkgreen" ,"darkgreen", "darkblue", "darkblue","darkred","darkgreen", "darkgreen", "darkgreen", "darkred", "darkred", "darkblue"), labelsize = 5, repel = TRUE, arrowsize = 1) + 
  theme(panel.grid = element_blank(),text = element_text(size = 18), axis.title = element_text(size = 20),
        axis.text = element_text(size = 20))


#PCA_clima___LAGOMERA----

pca.LAG<- read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/graficos corr-PCA/PCA_LIMPIOS CON ORLSS TABLA/PCA.lag.csv")
corr.LAG.cre <- as.data.frame(lapply(pca.LAG, as.numeric))
str(corr.LAG.cre)
corr_matrix <- cor(corr.LAG.cre)
head(corr_matrix)
ggcorrplot(corr_matrix)
p_value_matrix<- corr.test(corr.LAG.cre)$p #clacular matriz con p-valores
print(p_value_matrix)# Mostrar la matriz de p-valores
data.pca <- princomp(corr_matrix)
summary(data.pca)
data.pca$loadings[, 1:2] #cargar/mostrar componentes principales 1 y 2
fviz_eig(data.pca, addlabels = TRUE)

grp <- c("LN","PI", "P.HY", "SPEI.summer", "NAO.autumn","NAO.spring","Tme.summer", "Tme.spring" , "EA.spring")

fviz_pca_var(data.pca, col.var = grp, palette = c("darkblue","darkgreen", "darkred" ,"darkblue", "darkred", "darkgreen","darkred","darkblue", "darkblue"), labelsize = 5, repel = TRUE, arrowsize = 1) + 
  theme(panel.grid = element_blank(),text = element_text(size = 18), axis.title = element_text(size = 20),
        axis.text = element_text(size = 20)) 

#PCA_clima___LAPALMA----

pca.LAP <- read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/graficos corr-PCA/PCA_LIMPIOS CON ORLSS TABLA/PCA.lap.csv")
corr.LAP.cre <- as.data.frame(lapply(pca.LAP, as.numeric))
str(corr.LAP.cre)
corr_matrix <- cor(corr.LAP.cre)
head(corr_matrix)
ggcorrplot(corr_matrix)
p_value_matrix<- corr.test(corr.LAP.cre)$p #clacular matriz con p-valores
print(p_value_matrix)# Mostrar la matriz de p-valores
data.pca <- princomp(corr_matrix)
summary(data.pca)
data.pca$loadings[, 1:2] #cargar/mostrar componentes principales 1 y 2
fviz_eig(data.pca, addlabels = TRUE)

grp <- c("LN","PI", "P.autumn","NAO.autumn","NAO.winter", "P.annual" , "EA.spring")

fviz_pca_var(data.pca, col.var = grp, palette = c("darkblue","darkgreen", "darkred" ,"darkblue", "darkblue", "darkred","darkgreen"), labelsize = 3, repel = TRUE, arrowsize = 1) + 
  theme(panel.grid = element_blank(),text = element_text(size = 18), axis.title = element_text(size = 20),
        axis.text = element_text(size = 20)) 


###5.CORRELACIONES----

library('corrr')
library(ggcorrplot)
library("FactoMineR")
library(factoextra)
library(psych)


#------CORRELACIONES TENERIFE----

vcre.clima.TNF <-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/Datos climaticosP/corr.tenerife3.limp.csv")
str(vcre.clima.TNF)
numeric.cre.clima.TNF <- as.data.frame(lapply(vcre.clima.TNF, as.numeric))
str(numeric.cre.clima.TNF)
corr.TNF<- cor(numeric.cre.clima.TNF)
head(corr.TNF)
ggcorrplot(corr.TNF)

#write.table(corr.TNF,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\corr.TNF.csv")
print(corr.TNF)

#p-values

p_value_matrix_TNF <- corr.test(numeric.cre.clima.TNF)$p #clacular matriz con p-valores
print(p_value_matrix_TNF)

#write.table(p_value_matrix_TNF,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\p.valores.TNF.csv")


#------CORRELACIONES LA GOMERA----
vcre.clima.LAG <-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/Datos climaticosP/corr.lagomera3.limp.csv")
str(vcre.clima.LAG)
numeric.cre.clima.LAG <- as.data.frame(lapply(vcre.clima.LAG, as.numeric))
str(numeric.cre.clima.LAG)
corr.LAG<- cor(numeric.cre.clima.LAG)
head(corr.LAG)
ggcorrplot(corr.LAG)
#write.table(corr.LAG,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\corr.LAG.csv")
print(corr.LAG)

#p-values

p_value_matrix_LAG <- corr.test(numeric.cre.clima.LAG)$p #clacular matriz con p-valores
print(p_value_matrix)

#write.table(p_value_matrix_LAG,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\p.valores.LAG.csv")

#------CORRELACIONES LA PALMA----

vcre.clima.LAP <-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/Datos climaticosP/corr.lapalma3.limp.csv")
str(vcre.clima.LAP)
numeric.cre.clima.LAP <- as.data.frame(lapply(vcre.clima.LAP, as.numeric))
str(numeric.cre.clima.LAP)
corr.LAP<- cor(numeric.cre.clima.LAP)
head(corr.LAP)
ggcorrplot(corr.LAP)
#write.table(corr.LAP,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\corr.LAP.csv")
print(corr.LAP)

#p-values

p_value_matrix_LAP <- corr.test(numeric.cre.clima.LAP)$p #clacular matriz con p-valores
print(p_value_matrix)

#write.table(p_value_matrix_LAP,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\p.valores.LAP.csv")


#------CORRELACIONES AZORES----
vcre.clima.AZO<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/Datos climaticosP/corr.azores3.limp.csv")
str(vcre.clima.AZO)
numeric.cre.clima.AZO <- as.data.frame(lapply(vcre.clima.AZO, as.numeric))
corr_AZO <- cor(numeric.cre.clima.AZO) #calculo de matriz correlaciones
head(corr_AZO)
ggcorrplot(corr_AZO)

#p-values

p_value_matrix.AZO <- corr.test(numeric.cre.clima.AZO)$p #clacular matriz con p-valores
print(p_value_matrix.AZO)

#write.table(corr_AZO,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\corr.AZO.csv")
#write.table(p_value_matrix.AZO,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\p.valores.AZO.csv")

#------CORRELACIONES MADEIRA----

vcre.clima.MAD <-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/Datos climaticosP/corr.madeira3.limp.csv")
str(vcre.clima.MAD)
numeric.cre.clima.MAD <- as.data.frame(lapply(vcre.clima.MAD, as.numeric))
corr.MAD<- cor(numeric.cre.clima.MAD)
head(corr.MAD)
ggcorrplot(corr.MAD)
#write.table(corr.MAD,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\corr.MAD.csv")

#p-values

p_value_matrix_MAD <- corr.test(numeric.cre.clima.MAD)$p #clacular matriz con p-valores
print(p_value_matrix)

#write.table(p_value_matrix_MAD,"C:\\Users\\Usuario\\Desktop\\ARTICULO RESPUESTA CLIMA\\correlaciones\\p.valores.mad.csv")


# PLOTS CORRELACIONES CRECIMIENTO CLIMA----

corr_plot0 <-read.table("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/graficos corr-PCA/corr_plot_2.txt")

corr_plot0$MADLN <- as.numeric(corr_plot0$MADLN)
corr_plot0$MADOF <- as.numeric(corr_plot0$MADOF)
corr_plot0$MADCA <- as.numeric(corr_plot0$MADCA)
corr_plot0$TERIA <- as.numeric(corr_plot0$TERIA)
corr_plot0$TERJB <- as.numeric(corr_plot0$TERJB)
corr_plot0$TERLA <- as.numeric(corr_plot0$TERLA)
corr_plot0$TNFLN <- as.numeric(corr_plot0$TNFLN)
corr_plot0$TNFIP <- as.numeric(corr_plot0$TNFIP)
corr_plot0$TNFPE <- as.numeric(corr_plot0$TNFPE)
corr_plot0$TNFPI <- as.numeric(corr_plot0$TNFPI)
corr_plot0$TNFPL <- as.numeric(corr_plot0$TNFPL)
corr_plot0$LAGLN <- as.numeric(corr_plot0$LAGLN)
corr_plot0$LAGPI <- as.numeric(corr_plot0$LAGPI)
corr_plot0$LAPLN <- as.numeric(corr_plot0$LAPLN)
corr_plot0$LAPPI <- as.numeric(corr_plot0$LAPPI)

str(corr_plot0)

library(heatmaply)

heatmaply(corr_plot0,
          dendrogram = "none",
          margins = c(60,100,40,20),
          grid_color = "white",
          grid_width = 0.00001,
          fontsize_row = 15, fontsize_col = 15,
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient2(
            low = "blue", 
            high = "red", 
            midpoint = 0,
            limits = c(-0.5,0.5)
          ))



###6. Calcular Pointer Years con Point.Res----

#detr_crono <- detrend(crono, method = "Spline", nyrs = 30)

pyc <- pointer.norm(crono.MADLN , period = NULL, window = 13, method.thresh = "Cropper", C.thresh = 0.50, series.thresh = 50, make.plot = TRUE)
pyn <- pointer.norm(crono.TNFPL, period = NULL, window = 13, method.thresh = "Neuwirth", N.thresh = c(1, 1.28, 1.645), series.thresh = 50, make.plot = TRUE )
#event.plot(pyc, period = c(1940, 2010), x.tick.major = 10, x.tick.minor = 5)

cat("Pointer Years:", pyc$yr, "\n")

print(pyc)

rgc <- pointer.rgc(spline.crn.MADLN, period = NULL, nb.yrs = 4, rgc.thresh.pos = 50, rgc.thresh.neg = 50, series.thresh = 50, make.plot = TRUE)

it <- interval.trend(spline.crn.MADLN, period = NULL, trend.thresh = 0, IT.thresh = 95, make.plot = FALSE)
it<-interval.trend(spline.crn.AZLA, period = NULL, trend.thresh = 0, IT.thresh = 95, make.plot = TRUE)

detr_s033 <- detrend(s033, method = "Spline", nyrs = 30)
pyn <- pointer.norm(spline.crn.AZLA, method = "Neuwirth", make.plot = TRUE)

event.plot(pyn, period = c(1970, 2010), x.tick.major = 10, x.tick.minor = 5)

###7. INDICES RESISTENCIA/RECUPERACIÓN/RESILIENCIA___MACARONESIA (CON POINTER YEARS)----

resist_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/resistencia.macaronesia.PY.csv")
recov_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/recuperación.macaronesia.PY.csv")
resil_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/resiliencia.macaronesia.PY.csv")

resist <- as.data.frame(lapply(resist_chr, as.numeric))
recov <- as.data.frame(lapply(recov_chr, as.numeric))
resil <- as.data.frame(lapply(resil_chr, as.numeric))

str(resist)

boxplot(resist) + title("Resistance") 
boxplot(recov) + title("Recovery")  
boxplot(resil) + title("Resilience") 


###----INDICES RESISTENCIA/RECUPERACIÓN/RESILIENCIA___POR ESPECIE (CON POINTER YEARS)----

recov_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/INDICES POINTER YEARS/recuperacion.macaronesia.PY_especies.csv")
resist_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/INDICES POINTER YEARS/resistencia.macaronesia.PY_especies.csv")
resil_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/INDICES POINTER YEARS/resiliencia.macaronesia.PY_especies.csv")

resist <- as.data.frame(lapply(resist_chr, as.numeric))
recov <- as.data.frame(lapply(recov_chr, as.numeric))
resil <- as.data.frame(lapply(resil_chr, as.numeric))

str(resist)

nombres_especies<- c("Laurus azorica", "Juniperus brevifolia", "Ilex azorica", "Laurus novocanariensis", "Ocotea arborea", "Clethra arborea", "Persea indica", "Ilex perado", "Picconia excelsa", "Prunus lusitanica")
colores<- c("#24693D","#4A934E", "#7ABE6F" ,"#E6E600", "#F4D166", "#AD9024","#FEBEB1","#D0D0D0", "#858F93", "#515A65")


boxplot(resist, col = colores, main= "Recovery", cex.main=2, cex.axis=2) + 
  legend("topleft",legend = nombres_especies, fill = colores, cex = 1.2, box.lty = 0)


#par(mfcol = c(1, 3))
par(mfcol = c(3, 1))

boxplot(resist, col = colores, main= "Resistance", cex.main=2, cex.axis=2) 
boxplot(recov, col = colores, main= "Recovery", cex.main=2, cex.axis=2) 
boxplot(resil, col = colores, main= "Resilience", cex.main=2, cex.axis=2)

###----INDICES RESISTENCIA/RECUPERACIÓN/RESILIENCIA___POR ISLA (CON POINTER YEARS)----

recov_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/INDICES POINTER YEARS/recuperacion.macaronesia.PY_isla.csv")
resist_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/INDICES POINTER YEARS/resistencia.macaronesia.PY_isla.csv")
resil_chr<-read.csv("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/indices resiliencia/INDICES POINTER YEARS/resiliencia.macaronesia.PY_isla.csv")

resist <- as.data.frame(lapply(resist_chr, as.numeric))
recov <- as.data.frame(lapply(recov_chr, as.numeric))
resil <- as.data.frame(lapply(resil_chr, as.numeric))

str(resist)

nombres_islas<- c("Terceira", "Madeira","La Palma", "Tenerife", "La Gomera")
colores<- c("#499851","#E6E600","#D0D0D0","#F4D166","#FEBEB1")

#boxplot(resist, col = terrain.colors(ncol(resist))) + title("Resistance", cex.main=1.5) 

boxplot(resist, col = colores, main= "Recovery", cex.main=2, cex.axis=2) + 
  legend("topleft",legend = nombres_islas, fill = colores, cex = 1.2, box.lty = 0)

par(mfrow = c(1, 3))

boxplot(resist, col = colores, main= "Resistance", cex.main=2, cex.axis=2) 
boxplot(recov, col = colores, main= "Recovery", cex.main=2, cex.axis=2) 
boxplot(resil, col = colores, main= "Resilience", cex.main=2, cex.axis=2) 

boxplot(resil, col = terrain.colors(ncol(resist)), main= "Resilience", cex.main=2, cex.axis=2)

par(mfcol = c(3, 1))
boxplot(resist, col = terrain.colors(ncol(resist))) + title("Resistance") 
boxplot(recov, col = terrain.colors(ncol(resist))) + title("Recovery")  
boxplot(resil, col = terrain.colors(ncol(resist))) + title("Resilience")


##8. Models----


library(ggplot2)
library(lme4)
library(nortest)


#Recovery data

glmm.recov <- read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/modelo.recov.mac.csv")
str(glmm.recov)

plot(glmm.recov$recov)

glmm.recov$recov <- as.numeric(glmm.recov$recov)
glmm.recov$id <- as.factor(glmm.recov$id)
glmm.recov$species <- as.factor(glmm.recov$species)
glmm.recov$island <- as.factor(glmm.recov$island)
glmm.recov$year <- as.factor(glmm.recov$year)
glmm.recov$sqrt.recov <- sqrt(glmm.recov$recov)


#Resistance data

glmm.resist <- read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/modelo.resist.mac.csv")
str(glmm.resist)

plot(glmm.resist$resist)

glmm.resist$resist <- as.numeric(glmm.resist$resist)
glmm.resist$id <- as.factor(glmm.resist$id)
glmm.resist $species <- as.factor(glmm.resist$species)
glmm.resist$island <- as.factor(glmm.resist$island)
glmm.resist$year <- as.factor(glmm.resist$year)
glmm.resist$sqrt.resist <- sqrt(glmm.resist$resist)

#Resilience data

glmm.resil <- read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/modelo.resil.mac.csv")
str(glmm.resil)

plot(glmm.resil$resil)

glmm.resil$resil <- as.numeric(glmm.resil$resil)
glmm.resil$id <- as.factor(glmm.resil$id)
glmm.resil $species <- as.factor(glmm.resil$species)
glmm.resil$island <- as.factor(glmm.resil$island)
glmm.resil$year <- as.factor(glmm.resil$year)
glmm.resil$sqrt.resil <- sqrt(glmm.resil$resil)


hist(glmm.recov$recov)
hist(glmm.resist$resist)
hist(glmm.resil$resil)

hist(glmm.recov$sqrt.recov)
hist(glmm.resist$sqrt.resist)
hist(glmm.resil$sqrt.resil)

#AIC


#AIC(glmer(resist ~ island * species + (1|id) + (1|year), family = Gamma(link = "inverse"), glmm.resist))
#AIC(glmer(sqrt.resist ~ island * species + (1|id) + (1|year), family = Gamma(link = "inverse"), glmm.resist))
#AIC(glmer(sqrt.resist ~ island + (1|species) + (1|id) + (1|year), family = Gamma(link = "inverse"), glmm.resist))

AIC(glmer(resist ~ island + (1|species) + (1|id) , family = Gamma(link = "inverse"), glmm.resist)) 

#AIC(glmer(sqrt.recov ~ island + species + (1|id) + (1|year), family = Gamma(link = "inverse"), glmm.recov))
#AIC(glmer(sqrt.recov ~ island + (1|id) + (1|year), family = Gamma(link = "inverse"), glmm.recov))
#AIC(glmer(sqrt.recov ~ species + (1|id) + (1|year), family = Gamma(link = "inverse"), glmm.recov))

##RECOVERY----


model.recov.is <- glmer(sqrt.recov ~ island + (1|species) + (1|id) + (1|year),
                        family = Gamma(link = "inverse"), 
                        control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                        glmm.recov)


summary(model.recov.is)


library(emmeans)
library(multcompView)
library(multcomp)

emmeans_result <- emmeans(model.recov.is, ~ island)
cld_result <- multcomp::cld(emmeans_result, Letters = letters, adjust = "tukey")
print(cld_result)
cld_df <- as.data.frame(cld_result)
print(cld_df)
plot(cld_result)

model.recov.sp<-glmer(sqrt.recov ~ species + (1|island) + (1|id) + (1|year), 
                      family = Gamma(link = "inverse"), 
                      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                      glmm.recov)

boxplot(glmm.recov$recov ~ glmm.recov$species, las = 2 )

summary(model.recov.is)

##RESISTANCE----

model.resist.is<-glmer(sqrt.resist ~ island + (1|species) + (1|id)+ (1|year), 
                       control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                       family = Gamma, glmm.resist)

summary(model.resist.is)

emmeans_result <- emmeans(model.resist.is, ~ island)
cld_result <- multcomp::cld(emmeans_result, Letters = letters, adjust = "tukey")
print(cld_result)
cld_df <- as.data.frame(cld_result)
print(cld_df)
plot(cld_result)

resist<-posthoc (model.resist.is, EffectIndices = NULL, EffectLabels = NULL,
                 EffectsMatrix = NULL, ParBootstrap = FALSE, Nboots = 999,
                 SignificanceLevel = 0.05, UpperCase = FALSE,
                 RankLabels = TRUE, WaldApproximation = FALSE,
                 CalcClusters = FALSE, QUIET = TRUE, PlotAdj = FALSE,
                 digits = 4, padjust = NULL, Scale = 1.0, Location = 0.0,
                 isBinomialModel = FALSE, BackTransform = TRUE)

print(resist)

model.resist.sp<-glmer(resist ~ species + (1|island) + (1|id) + (1|year), 
                       family = Gamma, 
                       control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                       glmm.resist)

boxplot(glmm.resist$resist ~ glmm.resist$island, las = 2 )

summary(model.resist.sp)
sjPlot::tab_model(model.resist.is)
plot(glmm.resist$resist ,glmm.resist$sqrt.resist )
plot(glmm.resist$resist)


##RESILIENCE----


boxplot(glmm.resil$resil ~ glmm.resil$island, las = 2 )

model.resil.is<-glmer(resil ~ island + (1|species) + (1|id) + (1|year) ,
                      family = Gamma(link = "inverse"),
                      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                      glmm.resil)


emmeans_result <- emmeans(model.resil.is, ~ island)
cld_result <- multcomp::cld(emmeans_result, Letters = letters, adjust = "tukey")
print(cld_result)
cld_df <- as.data.frame(cld_result)
print(cld_df)
plot(cld_result)

resil<-posthoc (model.resil.is, EffectIndices = NULL, EffectLabels = NULL,
                EffectsMatrix = NULL, ParBootstrap = FALSE, Nboots = 999,
                SignificanceLevel = 0.05, UpperCase = FALSE,
                RankLabels = TRUE, WaldApproximation = FALSE,
                CalcClusters = FALSE, QUIET = TRUE, PlotAdj = FALSE,
                digits = 4, padjust = NULL, Scale = 1.0, Location = 0.0,
                isBinomialModel = FALSE, BackTransform = TRUE)

print(resil)

summary(model.resil.is)

model.resil.sp<-glmer(sqrt.resil ~ species + (1|island) + (1|id) + (1|year) ,
                      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                      family = Gamma(link = "inverse"), 
                      glmm.resil)

summary(model.resil.is)

sjPlot::tab_model(model.resil.is)

#MODEL- LAURUS----

#RECOVERY LAURUS----

recov_laurus<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/PI - LN/modelo.recov.mac-laurus.csv")
str(recov_laurus)

recov_laurus$recov <- as.numeric(recov_laurus$recov)
recov_laurus$id <- as.factor(recov_laurus$id)
recov_laurus$species <- as.factor(recov_laurus$species)
recov_laurus$island <- as.factor(recov_laurus$island)
recov_laurus$year <- as.factor(recov_laurus$year)
recov_laurus$sqrt <- sqrt(recov_laurus$recov)


model.recov.laurus<-glmer(sqrt ~ island + (1|id) + (1|year) ,
                          family = Gamma(link = "inverse"),
                          control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                          recov_laurus)


summary(model.recov.laurus)

pairs = emmeans(model.recov.laurus, ~ island)

pairs(pairs)
cld(pairs, Letters=letters)



#RESISTANCE LAURUS----

resist_laurus<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/PI - LN/modelo.resist.mac-laurus.csv")
str(resist_laurus)

resist_laurus$resist <- as.numeric(resist_laurus$resist)
resist_laurus$id <- as.factor(resist_laurus$id)
resist_laurus$species <- as.factor(resist_laurus$species)
resist_laurus$island <- as.factor(resist_laurus$island)
resist_laurus$year <- as.factor(resist_laurus$year)
resist_laurus$sqrt <- sqrt(resist_laurus$resist)


model.resist.laurus<-glmer(sqrt ~ island + (1|id) + (1|year) ,
                           family = Gamma(link = "inverse"),
                           control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                           resist_laurus)

summary(model.resist.laurus)

pairs = emmeans(model.resist.laurus, ~ island)

pairs(pairs)
cld(pairs, Letters=letters)

#RESILIENCE LAURUS----

resil_laurus<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/PI - LN/modelo.resil.mac-laurus.csv")
str(resil_laurus)

resil_laurus$resil <- as.numeric(resil_laurus$resil)
resil_laurus$id <- as.factor(resil_laurus$id)
resil_laurus$species <- as.factor(resil_laurus$species)
resil_laurus$island <- as.factor(resil_laurus$island)
resil_laurus$year <- as.factor(resil_laurus$year)
resil_laurus$sqrt <- sqrt(resil_laurus$resil)


model.resil.laurus<-glmer(sqrt ~ island + (1|id) + (1|year) ,
                          family = Gamma(link = "inverse"),
                          control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                          resil_laurus)

summary(model.resil.laurus)

pairs = emmeans(model.resil.laurus, ~ island)

pairs(pairs)
cld(pairs, Letters=letters)



#MODEL- PERSEA----

#RECOVERY PERSEA----

recov_persea<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/PI - LN/modelo.recov.mac-PI.csv")
str(recov_persea)

recov_persea$recov <- as.numeric(recov_persea$recov)
recov_persea$id <- as.factor(recov_persea$id)
recov_persea$species <- as.factor(recov_persea$species)
recov_persea$island <- as.factor(recov_persea$island)
recov_persea$year <- as.factor(recov_persea$year)
recov_persea$sqrt <- sqrt(recov_persea$recov)


model.recov.persea<-glmer(sqrt ~ island + (1|id) + (1|year) ,
                          family = Gamma(link = "inverse"),
                          control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                          recov_persea)


summary(model.recov.persea)



pairs = emmeans(model.recov.persea, ~ island)

pairs(pairs)
cld(pairs, Letters=letters)



#RESISTANCE PERSEA----

resist_persea<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/PI - LN/modelo.resist.mac-PI.csv")
str(resist_persea)

resist_persea$resist <- as.numeric(resist_persea$resist)
resist_persea$id <- as.factor(resist_persea$id)
resist_persea$species <- as.factor(resist_persea$species)
resist_persea$island <- as.factor(resist_persea$island)
resist_persea$year <- as.factor(resist_persea$year)
resist_persea$sqrt <- sqrt(resist_persea$resist)


model.resist.persea<-glmer(sqrt ~ island + (1|id) + (1|year) ,
                           family = Gamma(link = "inverse"),
                           control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                           resist_persea)

summary(model.resist.persea)

pairs = emmeans(model.resist.persea, ~ island)

pairs(pairs)
cld(pairs, Letters=letters)

#RESILIENCE PERSEA----

resil_persea<-read.csv2("C:/Users/Usuario/Desktop/ARTICULO RESPUESTA CLIMA/modelos mixtos/PI - LN/modelo.resil.mac-PI.csv")
str(resil_persea)

resil_persea$resil <- as.numeric(resil_persea$resil)
resil_persea$id <- as.factor(resil_persea$id)
resil_persea$species <- as.factor(resil_persea$species)
resil_persea$island <- as.factor(resil_persea$island)
resil_persea$year <- as.factor(resil_persea$year)
resil_persea$sqrt <- sqrt(resil_persea$resil)


model.resil.persea<-glmer(sqrt ~ island + (1|id) + (1|year) ,
                          family = Gamma(link = "inverse"),
                          control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)),
                          resil_persea)

summary(model.resil.persea)

pairs = emmeans(model.resil.persea, ~ island)

pairs(pairs)
cld(pairs, Letters=letters)


###_4_Assumpitions of the models----


#Residuals----

residuos.recov<-residuals(model.recov.is) # estandarizados? 
plot(residuos.recov) 
hist(residuos.recov, xlab = "Residuals", main = "")


residuos.resist<-residuals(model.resist.is) # estandarizados? 
plot(residuos.resist) 
hist(residuos.resist, xlab = "Residuals", main = "")


residuos.resil<-residuals(model.resil.is) # estandarizados? 
plot(residuos.resil) 
hist(residuos.resil, xlab = "Residuals", main = "")


#Homocedasticity---- 

island.recov <- glmm.recov$island
plot(island.recov, residuos.recov, xlab = "site", ylab = "Residuals")


island.resist <- glmm.resist$island
plot(island.resist, residuos.resist, xlab = "site", ylab = "Residuals")


island.resil <- glmm.resil$island
plot(island.resil, residuos.resil, xlab = "site", ylab = "Residuals")


#Residuals vs predicted----

fit.recov <- fitted(model.recov.is)
plot(residuos.recov ~ fit.recov, ylab="Residuals", xlab="Fitted")

fit.resist <- fitted(model.resist.is)
plot(residuos.resist ~ fit.resist, ylab="Residuals", xlab="Fitted")

fit.resil <- fitted(model.resil.is)
plot(residuos.resil ~ fit.resil, ylab="Residuals", xlab="Fitted")


#Plots----

par(mfcol=c(2,2))

hist(residuos.recov, xlab = "Residuals", main = "")
plot(island.recov, residuos.recov, xlab = "site", ylab = "Residuals")
plot(residuos.recov ~ fit.recov, ylab="Residuals", xlab="Fitted")

plot(model.recov.is)


par(mfcol=c(2,2))

hist(residuos.resist, xlab = "Residuals", main = "")
plot(island.resist, residuos.resist, xlab = "site", ylab = "Residuals")
plot(residuos.resist ~ fit.resist, ylab="Residuals", xlab="Fitted")

plot(model.resist.is)


par(mfcol=c(2,2))

hist(residuos.resil, xlab = "Residuals", main = "")
plot(island.resil, residuos.resil, xlab = "site", ylab = "Residuals")
plot(residuos.resil ~ fit.resil, ylab="Residuals", xlab="Fitted")

plot(model.resil.is)