Práctico 7. Selección de Modelos

library(arm)
library(car)
library(GGally)
library(glmnet)
library(leaps)
library(MASS) 
library(MuMIn)
library(performance)

dat <- read.table("s_poly.txt", header = TRUE)

# exploración de los datos
layout(matrix(1:8,2,4))
plot(dat$azucar); plot(dat$largo.l.tot); plot(dat$largo.tubo) 
plot(dat$caliz.sup); plot(dat$caliz.med); plot(dat$caliz.inf) 
plot(dat$labio.sup); plot(dat$labio.inf)
layout(1)

# exploración bivariada de los datos
pairs(dat)

# Detección de la colinealidad
# A) Matrices de correlación
CORR <- cor(dat[, c(
        "largo.l.tot", "largo.tubo", "caliz.sup",
        "caliz.med", "caliz.inf", "labio.sup", "labio.inf")], 
        use = "complete.obs")

# B) scatterplot matrix
ggpairs(dat, aes(alpha = 0.4), 
        upper = list(continuous = wrap('cor', size = 3)))

# cuando aumenta el número de variables, puede ser útil
# no visualizar los datos individuales
ggcorr(dat, method = c("pairwise", "pearson"))

# C) Factores de inflación de la varianza
fit1 <- lm(azucar ~ largo.l.tot + largo.tubo + caliz.sup + caliz.med +
        caliz.inf + labio.sup + labio.inf, data = dat)

vif(fit1)   
check_collinearity(fit1)
ckeck_model(fit1)

# Para debatir: descartamos el cáliz inferior y construimos un nuevo modelo
fit2 <- lm(azucar ~ largo.l.tot + largo.tubo + caliz.sup + caliz.med +
        labio.sup + labio.inf, data = dat)

vif(fit2)
check_model(fit)

summary(fit2)
anova (fit2)


# LIBRARY leaps
# Selección de modelos. Best subset selection.
# Construyendo todos los modelos posibles
# (pueden utilizarse forward o backward cambiando  el argumento method)

BS.fit <- regsubsets(azucar ~ largo.l.tot + largo.tubo + caliz.sup +
        caliz.med + labio.sup + labio.inf, data = dat, method = "exhaustive")
BS.summary <- summary(BS.fit)
BS.summary

plot(BS.summary$cp, xlab = "número de variables", ylab = "Cp", type ="l")
plot(BS.summary$bic, xlab = "número de variables", ylab = "BIC", type ="l")

plot(BS.fit, scale = "Cp")
plot(BS.fit, scale = "bic")

coef(BS.fit, 3)
coef(BS.fit, 1)

# LIBRARY MASS
# stepAIC realiza una búsqueda automática del mejor modelo 
stepAIC(fit2, direction = "both") # para AIC
n <- nrow(dat)
stepAIC(fit2, direction = "both", k = log(n)) # para BIC

# LIBRARY Arm; MuMIn
# Model averaging
global <- lm(azucar ~ largo.l.tot + largo.tubo + caliz.sup +
        caliz.med + labio.sup + labio.inf, data = dat, na.action = na.fail)
std.model <- standardize(global, standardize.y = FALSE)
set <- dredge(std.model)
set

top.mod <- get.models(set, subset = delta < 2)
model.sel(top.mod)

AVG <- model.avg(top.mod)
summary(AVG)
sw(AVG)

# LIBRARY glmnet
# Selección de Modelos. Shrinkage - Lasso

x <- model.matrix(azucar ~ largo.l.tot + largo.tubo + caliz.sup +
        caliz.med + labio.sup + labio.inf, data = dat)[, -1]
y <- dat$azucar

lasso.mod <- glmnet(x, y, alpha = 1)
plot(lasso.mod)
coef(lasso.mod, s = 0.1) # ejemplo

# Cross validation.
cv.out <- cv.glmnet(x, y, alpha = 1, nfolds = 10)
plot(cv.out)

cv.out$lambda.min
coef(cv.out, s = "lambda.min")