# Capítulo 9. Extensiones de la regresión

#  9.1. Modelos polinómicos

library(ADER)
data(dry)
str(dry)
     
glm.dry <- glm(richness~ndvi, data=dry,
  family = poisson)

dry$ndvi_poly2 <- dry$ndvi^2
glm.poly2 <- glm(richness~ndvi + ndvi_poly2, data=dry,
  family = poisson)

glm.poly2 <- glm(richness~ndvi + I(ndvi^2), dry,
  family = poisson)

glm.poly2ortho <- glm(richness~poly(ndvi, degree=2), dry,
  family = poisson)

dry$ndvi_poly2ortho <- (dry$ndvi - mean(dry$ndvi))^2

cor(dry$ndvi, dry$ndvi_poly2)
cor(dry$ndvi, dry$ndvi_poly2ortho)

library(car)
vif(glm.poly2)
anova(glm.poly2, test="Chi")
anova(glm.poly2ortho, test="Chi")

glm.poly3ortho <- glm(richness~poly(ndvi, degree = 3),
  dry, family = poisson)
anova(glm.dry, glm.poly2ortho, glm.poly3ortho, test="Chi")

plot(richness~ndvi, data=dry,pch=16, col="grey",
  xlab="Productividad media (NDVI)", ylab = "Riqueza de árboles")

xv <- seq(0.4, 0.9, 0.001)
yv2 <- predict(glm.poly2ortho, list(ndvi = xv), "response")
yv3 <- predict(glm.poly3ortho, list(ndvi = xv), "response")
lines(xv, yv2, lwd=2)
lines(xv, yv3, lwd=2, lty=3)

library(DHARMa)
simulationOutput <- simulateResiduals(glm.poly2ortho, n = 1000)
plot(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput)

library(MASS)
nb.poly2ortho <- glm.nb(richness~poly(ndvi, 2), dry)
simulationOutput <- simulateResiduals(nb.poly2ortho, n = 1000)
plot(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput)

library(modEvA)
Dsquared(nb.poly2ortho)

#  9.2. Modelos aditivos generalizados (GAM)
#  9.2.3. Funciones para ajustar un GAM en R
#  A) El paquete gam

library(gam)
gam1 <- gam(richness~s(ndvi), data=dry, family=poisson)
summary(gam1)

par(mfcol=c(1,2))
y_hat <- predict(gam1, se.fit=TRUE)
ord <- order(dry$ndvi)
plot(richness~ndvi, data=dry, las=1,
  xlab="Productividad media (NDVI)", ylab="Riqueza")
lines(dry$ndvi[ord], exp(y_hat$fit)[ord], lwd=2)
lines(dry$ndvi[ord], exp(y_hat$fit + 2*y_hat$se.fit)[ord], lty=2)
lines(dry$ndvi[ord], exp(y_hat$fit - 2*y_hat$se.fit)[ord], lty=2)

plot(gam1, xlab="Productividad media (NDVI)",
  ylab="s(NDVI)", se=TRUE)

library(modEvA)
Dsquared(gam1)
gam_1n <- gam(richness~s(ndvi, df=1), data=dry, family=poisson)
gam_2n <- gam(richness~s(ndvi, df=2), data=dry, family=poisson)
gam_4n <- gam(richness~s(ndvi, df=4), data=dry, family=poisson)
gam_8n <- gam(richness~s(ndvi, df=8), data=dry, family=poisson)

library(gam)
gam2 <- gam(richness~lo(ndvi), data=dry, family=poisson)
summary(gam2)
coefficients(gam2)

library(gam)
gam2 <- gam(richness~lo(ndvi), data=dry, family=poisson)
Dsquared(gam2)
plot(gam2, se=T, xlab="Productividad media (NDVI)", ylab="l(NDVI)")

gam_1s <- gam(richness~lo(ndvi, span=0.1), data=dry, family=poisson)
gam_2s <- gam(richness~lo(ndvi, span=0.2), data=dry, family=poisson)
gam_5s <- gam(richness~lo(ndvi, span=0.5), data=dry, family=poisson)
gam_9s <- gam(richness~lo(ndvi, span=0.9), data=dry, family=poisson)

#  B) El paquete mgcv

library(mgcv)
mgcv.gam1 <- mgcv::gam(richness~s(ndvi), data=dry, family=poisson)
summary(mgcv.gam1)

par(mfcol=c(1,2))
y_hat <- predict(mgcv.gam1, se.fit=TRUE)
ord <- order(dry$ndvi)
plot(richness~ndvi, data=dry, las=1,
  xlab="Productividad media (NDVI)", ylab="Riqueza")
lines(dry$ndvi[ord], exp(y_hat$fit)[ord], lwd=2)
lines(dry$ndvi[ord], exp(y_hat$fit + 2*y_hat$se.fit)[ord], lty=2)
lines(dry$ndvi[ord], exp(y_hat$fit - 2*y_hat$se.fit)[ord], lty=2)

plot(mgcv.gam1, xlab="Productividad media (NDVI)",
  ylab="s(NDVI)")

mgcv.gam4k <- mgcv::gam(richness~s(ndvi, k=4, fx=TRUE), data=dry,
  family=poisson)
plot(mgcv.gam4k)

#  9.2.4. Evaluación del modelo

par(mfcol=c(2,2))
gam.check(mgcv.gam1, pch=19)
simulationOutput <- simulateResiduals(mgcv.gam1)
plot(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput)

mgcv.gam2 <- mgcv::gam(richness~s(ndvi), data=dry, family=nb())
summary(mgcv.gam2)
par(mfcol=c(2,2))
gam.check(mgcv.gam2, pch=19)

library(aods3)
gof(mgcv.gam1)
gof(mgcv.gam2)

#  9.3. Árboles de regresión y de clasificación (CART)

#  9.3.1. El paquete rpart

library(ADER)
data(defo)
str(defo)

library(rpart)
set.seed(2)
rpart.def <- rpart(Defoliacion~Especie+Densidad_pinos+Elevacion+
  Orientacion+Insolacion+Potencial_hidrico, defo)
rpart.def

library(rpart.plot)
rpart.plot(rpart.def, type=3)

par(mfcol=c(1,2))
plot(predict(rpart.def), residuals(rpart.def),
  xlab="Valores esperados", ylab="Residuos")
abline(h=0)
qqnorm(residuals(rpart.def), xlab="Cuantiles teóricos",
  ylab="Cuantiles teóricos")
qqline(residuals(rpart.def))

printcp(rpart.def)

par(mfrow=c(1,2))
rsq.rpart(rpart.def)
minXE <- min(rpart.def$cptable[,"xerror"])
WminXE <- which.min(rpart.def$cptable[,"xerror"])
min.xstd<- rpart.def$cptable[WminXE,"xstd"]
abline(h=minXE +min.xstd, lty=3)

prune.rpart.def <- prune.rpart(rpart.def, cp=0.0759)
rpart.plot(prune.rpart.def, type=3, cex=1.5)

defo$categ <- cut(defo$Defoliacion, breaks=c(0,25,50,100),
  labels=c("Baja", "Media", "Alta"))

rpart.categ <- rpart(categ~Especie+Densidad_pinos+Elevacion+
  Orientacion+Insolacion+Potencial_hidrico, defo)
rpart.plot(rpart.categ, type=3)

#  9.3.2. Alternativas a los árboles de regresión y clasificación

library(party)
plot(ctree(Defoliacion~Especie+Densidad_pinos+Elevacion, data=defo))

library(party)
plot(mob(Defoliacion~ Elevacion | Especie+Densidad_pinos+
  Orientacion+Insolacion+Potencial_hidrico, data=defo))

library(randomForest)
set.seed(3)
bgg.def <- randomForest (Defoliacion~Especie+Densidad_pinos+
  Elevacion+Orientacion+Insolacion+Potencial_hidrico, data=defo,
  mtry=6, importance=TRUE)
importance(bgg.def)

library(randomForest)
set.seed(3)
rf.def <- randomForest (Defoliacion~Especie+Densidad_pinos+
  Elevacion+Orientacion+Insolacion+Potencial_hidrico, data=defo,
  importance=TRUE)
varImpPlot(rf.def)

library(gbm)
set.seed(3)
gbm.def <- gbm(Defoliacion~Especie+Densidad_pinos+Elevacion+
  Orientacion+Insolacion+Potencial_hidrico, data=defo,
  distribution = "gaussian", n.trees =5000, interaction.depth =4)

par(mar=c(5,10,4,2))
summary(gbm.def, las=1)