# Capítulo 4. Introducción a los modelos lineales: regresión

#  4.2. Regresión simple
#  4.2.2. Cómo ajustar un modelo lineal en R

library(ADER)
data(cladonia)
plot(cladrang~Humdepth, data=cladonia,
  ylab="Cobertura de Cladonia (%)",
  xlab="Grosor de la capa de humus del suelo (cm)")
  
lm.clad_H <- lm(cladrang~Humdepth, data=cladonia)

#  4.2.3. Contraste de hipótesis y variación explicada por el
#           modelo

anova(lm.clad_H)
sum((cladonia$cladrang - mean(cladonia$cladrang))^2)
sum(anova(lm.clad_H)$"Sum Sq")


R2 <- anova(lm.clad_H)$"Sum Sq"[1]/sum(anova(lm.clad_H)$"Sum Sq")
R2


#  Cuadro 4.1. Estimación por mínimos cuadrados

meanX <- mean(cladonia$Humdepth)
meanY <- mean(cladonia$cladrang)
n <- length(cladonia$cladrang)
var.xy <- sum((cladonia$cladrang - meanY)*
  (cladonia$Humdepth - meanX))/(n-1)
var.x2 <- sum((cladonia$Humdepth - meanX)^2)/(n-1)
(B1 <- var.xy/var.x2)
(B0 <- meanY-(B1*meanX))

coefficients(lm.clad_H)

#  4.2.4. Interpretación del modelo

summary(lm.clad_H)
plot(cladrang~Humdepth, data=cladonia,
  ylab="Cobertura de Cladonia (%)",
  xlab="Grosor de la capa de humus del suelo (cm)")

lm.clad_H <- lm(cladrang~Humdepth, data=cladonia)
abline(lm.clad_H, lwd=2)
plot(cladrang~Humdepth, data=cladonia,
  ylab="Cobertura de Cladonia (%)",
  xlab="Grosor de la capa de humus del suelo (cm)")

lm.clad_H <- lm(cladrang~Humdepth, data=cladonia)

xv <- seq(0, 4, 0.01)
predict.clad <- predict(lm.clad_H, newdata=list(Humdepth=xv),
  interval="confidence")
lines(xv, predict.clad[,1], lwd=2)
lines(xv, predict.clad[,2], lwd=2, lty=3)
lines(xv, predict.clad[,3], lwd=2, lty=3)


#  Cuadro 4.2. Sesgos en el R2 cuando hay muchas
#                   variables explicativas

y <- rnorm(1000)
df <- data.frame(matrix(rnorm(990*1000), ncol=990, nrow=1000))
results <- data.frame(R2=NULL, R2adj=NULL)
index <- c(1, seq(10, 990, 10))

for(i in 1:length(index)) {
   lmfit <- lm(y ~ ., df[1:index[i]])
   results2 <- data.frame(R2=summary(lmfit)$r.squared,
   R2adj=summary(lmfit)$adj.r.squared)
   results <- rbind(results, results2)
}

plot(results$R2~index, ylim=c(-0.1,1),
  ylab="R2 y R2 ajustado",
  xlab="N. variables en el modelo", pch=16)

points(index, results$R2adj, pch=4)
abline(h=0)
legend("topleft", pch=c(16,4), legend=c("R2", "R2 ajustado"))

#  4.2.5. Revisión de los supuestos del modelo

par(mfrow=c(2,2))
plot(lm.clad_H)

shapiro.test(residuals(lm.clad_H))

library(lmtest)
resettest(lm.clad_H)

library(car)
ncvTest(lm.clad_H)

library(car)

library(lmtest)
coeftest(lm.clad_H, vcov=hccm)


#  Cuadro 4.3. ¿Qué se puede hacer cuando no se
#                  cumplen los supuestos del modelo?

library(faraway)
data(gala)
par(mfrow=c(2,2))

lm.gala <- lm(Species~Area, gala)
lm.log <- lm(log(Species)~log(Area), gala)

plot(Species ~ Area, gala)
abline(lm.gala)

plot(log(Species)~log(Area), gala)
abline(lm.log)

plot(residuals(lm.gala)~fitted(lm.gala))
abline(h=0)

plot(residuals(lm.log)~fitted(lm.log))
abline(h=0)

#  4.3. Regresión múltiple
#  4.3.1. Ajuste del modelo

library(ADER)
data(cladonia)
variab <- c("cladrang", "N", "P", "Humdepth")

library(SciViews)
pairs(cladonia[, variab],
  upper.panel=panel_cor,
  lower.panel=panel.reg)
lm.clad <- lm(cladrang~N+P+Humdepth, cladonia)

#  4.3.2. Contraste de hipótesis y variación explicada por el
#            modelo

anova(lm.clad)

#   C) Especificar diferentes sumas de cuadrados en R

library(car)
Anova(lm.clad, type="III")

#  4.3.3. Interpretación y representación gráfica del modelo
#  ajustado

summary(lm.clad)
n <- 100
nuevos.datos <- data.frame(N=seq(min(cladonia$N)-5,
  max(cladonia$N)+5, le=n), P=rep(mean(cladonia$P), n),
  Humdepth=rep(mean(cladonia$Humdepth), n))

respuesta.N <- predict(lm.clad, newdata=nuevos.datos,
  type="response")

with(cladonia, plot(N, cladrang, main="",
  ylab="Cobertura de Cladonia (%)", xlab="Concentración de N"))
lines(nuevos.datos$N, respuesta.N)

library(effects)
plot(allEffects(lm.clad), ci.style="lines", row=1, col=3,
  cex.lab=1.1, cex.axis=1.1)

library(visreg)
par(mfrow=c(1,3))
visreg(lm.clad)

#  4.3.4. Coeficientes de regresión estandarizados

lm.clad.st <- lm(scale(cladrang)~scale(N)+scale(P)+
  scale(Humdepth),cladonia)
  
summary(lm.clad.st)

#  4.3.5. Multicolinealidad

variab <- c("N", "P", "K", "Ca", "Mg", "S", "Humdepth", "pH")
library(corrgram)
corrgram(cladonia[, variab],
  lower.panel=panel.pts,
  upper.panel=panel.conf,
  diag.panel=panel.density)

library(car)
lm.clad.comp <- lm(cladrang~., data=cladonia)
vif(lm.clad.comp)


#  Cuadro 4.4. Regresión parcial y partición de la
#                    varianza

m.abc <- lm(cladrang~N+P+Humdepth, data=cladonia)
(p.abc <- summary(m.abc)$adj.r.squared)

m.ab <- lm(cladrang~N+P, data=cladonia)
(p.ab <- summary(m.ab)$adj.r.squared)

m.bc <- lm(cladrang~Humdepth, data=cladonia)
(p.bc <- summary(m.bc)$adj.r.squared)
(p.a <- p.abc-p.bc)
(p.c <- p.abc - p.ab)
(p.b <- p.bc-p.c)

library(yhat)
regr(lm.clad)$Commonality_Data$CC

library(vegan)
(clad.vp <- varpart(cladonia$cladrang, ~N+P, ~Humdepth,
  data=cladonia))
plot(clad.vp, Xnames=c("N+P", "Humdepth"))

#  4.4. Regresión de tipo I y de tipo II

library(smatr)
data(leaflife)
lm.leaf <- lm(log10(longev)~log10(lma), data=leaflife)
ma.leaf <- ma(longev~lma, log='xy', data=leaflife)
sma.leaf <- sma(longev~lma, log='xy', data=leaflife)

plot(log10(longev)~log10(lma), leaflife,
  xlab="Log(Peso específico de hoja)",
  ylab="Log(Longevidad de hoja)")
abline(lm.leaf, lwd=2)
abline(ma.leaf, lwd=2, lty=3)
abline(sma.leaf, lwd=2, lty=2)

legend("bottomright", lty=c(1,3,2), lwd=rep(2,3),
legend=c("Modelo I", "Modelo II (MA)", "Modelo II (SMA)"))