Se intenta realizar un modelo que prediga el contenido de azúcar en el néctar en flores de Salvia polystachia, a partir de 7 variables morfológicas medidas. Los datos se encuentran en el archivo s_poly.txt.
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")