#################################################################
##                                                             ##
##            MODELOS GAM - MACROINVERTEBRADOS                 ##
##                Neus Sepulcre Vidal                          ##
##                                                             ##
#################################################################

gam <- read.table("C:/gam.txt", sep = "\t", header=T)
names(gam)
[1] "tramo"            "region"           "Elevación"        "Temperatura"      "Caudal"           "calado"           "velocidad"        "caladoxvelocidad"
[9] "Algas"            "Moss"             "Macrophyt"        "Bedrock"          "Sand"             "Gra1.4"           "Gra5.8"           "Gra9.12"         
[17] "Granulom."        "Leuctra"          "Rhyacophila"      "Elmis"            "Simuliidae"       "Ceratopogonidae"  "Hydraena"         "Protonemura"     
[25] "Ostracoda"        "Rhithrogena"      "Sericostoma"      "Ephemerella"      "Tinodes"          "Polycelis" 

attach(gam)
tramo <- cbind (tramo)
region <- cbind (region)
Elevación <- cbind (Elevación)
Temperatura <- cbind (Temperatura)
Caudal <- cbind (Caudal)
calado <- cbind (calado)
velocidad <- cbind (velocidad)
caladoxvelocidad <- cbind (caladoxvelocidad)
Algas <- cbind (Algas)
Moss <- cbind (Moss)
Macrophyt <- cbind (Macrophyt)
Bedrock <- cbind (Bedrock)
Sand <- cbind (Sand)
Gra1.4 <- cbind (Gra1.4)
Gra5.8 <- cbind (Gra5.8)
Gra9.12 <- cbind (Gra9.12)
Granulom. <- cbind (Granulom.)
Leuctra <- cbind (Leuctra)
detach (gam)



#############################################################
##                                                         ##
##      MODELOS GAM -                                      ##
##                                                         ##
#############################################################

Conjunto.Leuctra <-cbind (Leuctra,tramo,region,Elevación,Temperatura,Caudal,calado,velocidad,caladoxvelocidad,Algas,Moss,Macrophyt,Bedrock,Sand,Gra1.4,Gra5.8,Gra9.12,Granulom.)    # Junto las variables para crear la tabla de trabajo
Conjunto.Leuctra <- as.data.frame (Conjunto.Leuctra)    # Creación de tabla

# Observo las correlaciones de las variables del modelo
# Con curvas de suavizado, histogramas de las variables y coef de correlación
# Con esto observo la correlación entre variables. 0 > 0.05 las elimino del modelo

    panel.hist <- function(x, ...)
    {
       usr <- par("usr"); on.exit(par(usr))
       par(usr = c(usr[1:2], 0, 1.5) )
       h <- hist(x, plot = FALSE)
       breaks <- h$breaks; nB <- length(breaks)
       y <- h$counts; y <- y/max(y)
       rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
    }
    panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
    {
       usr <- par("usr"); on.exit(par(usr))
       par(usr = c(0, 1, 0, 1))
       r <- cor(x, y, use="complete")
       txt <- format(c(r, 0.123456789), digits=digits)[1]
       txt <- paste(prefix, txt, sep="")
       if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
       text(0.5, 0.5, txt, cex = cex.cor * r)
    }
pairs(Conjunto.Leuctra, lower.panel=panel.smooth,upper.panel=panel.cor,diag.panel=panel.hist)
dev.off()

# Antes de pasar a los GAM miro las variables lineales más significativas mediante GLM
# Selección Stepwise basado AIC
library(MASS)
modelo.1 <-glm(Leuctra~.,data=Conjunto.Leuctra)
modelo.1.step <- stepAIC(modelo.1,direction="both",trace=T)

# MODELOS ADITIVOS GENERALIZADOS (GAMs)
#--------------------------------------

# Modelos aditivos usando mgcv()
#---------------------------------

library(mgcv)
library(nlme)
library(gam)

# Los grados de libertad para la función de suavizado está elegido internamente por defecto. El paquete gam depende del usuario.

Leuctra.mgcv.1<-gam(Leuctra ~ s(Dadul) +s(ANCH) + s(PROFmed) + s(Sf), data=Conjunto.Leuctra)
summary(Leuctra.mgcv.1)

# Examinar el modelo mediante gráficos
par(mfrow=c(4,1))   # Poner en 4 el número de variables que te salgan del modelo gam
plot(Leuctra.mgcv.1,residuals=F,se=T,pch=".",all.terms=T,shade=TRUE,shade.col="#A6CEE3",seWithMean=TRUE)
dev.off()

# Diagnóstico del paquete mgcv
qq.gam(Leuctra.mgcv.1,rep=50,level=1,type=c("deviance","pearson","response"),pch=19,cex=.2)
gam.check (Leuctra.mgcv.1, old.style=FALSE,type=c("deviance","pearson","response"),rep=50, level=.9, rl.col=2, rep.col="gray80")
dev.off()

# con presencia/ausencia la distribución es binomial
binomial.mgcv.Leuctra <-gam(Leuctra ~ s(Dadul) +s(ANCH) + s(PROFmed) + s(Sf), family=binomial(link = "logit"), data=Conjunto.Leuctra)
summary(binomial.mgcv.Leuctra)

par(mfrow=c(4,1))
plot(binomial.mgcv.Leuctra, residuals=F,se=T,pch=".",all.terms=T,shade=TRUE,shade.col="#A6CEE3",seWithMean=TRUE)
dev.off()

vis.gam(Leuctra.mgcv.1,theta=35,phi=20,color="gray", plot.type="persp")
vis.gam(Leuctra.mgcv.1,theta=35,color="gray", plot.type="contour")
dev.off()

# Validación y comparación del modelo
#------------------------------------

# Generalized cross validation (GCV) criterion
Leuctra.mgcv.1$gcv.ubre

#Akaike information criterion (AIC)
Leuctra.mgcv.1$aic