Base de datos

En este documento se utlizará la base de datos mpg de la librería gamair, la cual contiene información sobre el consumo de combustible de una variedad de vehículos en USA, con sus respectivas caraterísticas mecánicas y físicas. La base de datos tiene 205 observaciones y 26 variables.

require(gamair)


getdata <- function(...)    ##Función para obtener datos de paquetes
{
    e <- new.env()        
    name <- data(..., envir = e)[1]
    e[[name]]
}

datos <- getdata("mpg")    ##Guardando los datos
attach(datos)              ##Código para acceder a los datos con el nombre

Objetivo

El objetivo de este documento será ajustar algunos modelos generalizados aditivos(GAMs), haciendo uso del paquete mgcv con la base de datos mpg, para las variables respuesta city.mpg y hw.mpg, las cuales representan el consumo de combustible en millas por galón en la ciudad y el consumo de combustible en millas por galón en carretera, respectivamente.

Análisis descriptivo

Como primer paso se presenta información sobre las variables respuesta.

library(ggplot2)
library(ggpubr)

##Construcción de los histogramas
g0 <- ggplot(datos, aes(x = city.mpg)) +
  geom_histogram(aes(y =..density..), color = 'black', fill = 'salmon', 
                 binwidth = 3) +
  labs(title = "Fuel consumption in town as miles per gallon.", x = "Miles per gallon ", y = "Densidad") + 
  theme_bw()
  
  
g1 <- ggplot(datos, aes(x = hw.mpg)) + 
  geom_histogram(aes(y =..density..), color = 'black', fill = 'aquamarine', binwidth = 3) + 
  labs(title = "Fuel consumption on highway as miles per US gallon.", x = "Miles per gallon ", y = "Densidad") + 
  theme_bw()

##Graficando los histogramas en una sola página
ggarrange(g0, g1, ncol = 2, nrow = 1)

Luego de tener el histograma de las variables, se intentará buscar la distribución que mejor se ajuste a los datos, de manera que el modelo se ajuste con la distribución de la familia exponencial más adecuada.

library(fitdistrplus)          ##Paquete para ajustar distribuciones a los datos
descdist(city.mpg, discrete = FALSE) ##Función para encontrar mejor ajuste

## summary statistics
## ------
## min:  13   max:  49 
## median:  24 
## mean:  25.21951 
## estimated sd:  6.542142 
## estimated skewness:  0.663704 
## estimated kurtosis:  3.578648
library(fitdistrplus)          ##Paquete para ajustar distribuciones a los datos
descdist(hw.mpg, discrete = FALSE) ##Función para encontrar mejor ajuste

## summary statistics
## ------
## min:  16   max:  54 
## median:  30 
## mean:  30.75122 
## estimated sd:  6.886443 
## estimated skewness:  0.5399972 
## estimated kurtosis:  3.44007

La curtosis y el cuadrado de la misma, muestran que ambas variables pueden distribuirse gamma, normal, lognormal o inclusive weibull, pero para efectos de este ejercicio se usará la distribución normal y alguno de los modelos que se plantearan se ajustará además con una distribución gamma.

city.norm <- fitdist(city.mpg, "norm") ##Ajustando los datos a una distribución normal
plot(city.norm) ##Chequeo gráfico del ajuste

hw.norm <- fitdist(city.mpg, "norm")##Ajustando los datos a una distribución normal
plot(hw.norm) ##Chequeo gráfico del ajuste

Las figuras anteriores muestran que los datos siguen de cierta manera una distribución normal.

Ahora se selecionan algunas variables del conjunto de datos que serán usadas para entrenar el modelo, estás variables son:
weight: el peso del vehículo en libras.

hp: caballos de fuerza del vehículo.

eng.cc: capacidad del motor en pulgadas cúbicas.

##Creando el conjunto de datos para cada variable independiente

myvars.city <- c("weight", "eng.cc", "hp", "city.mpg")
data.city <- datos[myvars.city]

myvars.hw <- c("weight", "eng.cc", "hp", "hw.mpg")
data.hw <- datos[myvars.hw]

##Gráfico de correlaciones con ajuste LOESS
library(psych)

pairs.panels(data.city, 
             method = "pearson", # Método para el cálculo de las correlaciones
             hist.col = "aquamarine",
             density = TRUE,  # Mostrando las densidades
             ellipses = FALSE #Ignorando las elipses de correlación
             )

##Gráfico de correlaciones con ajuste LOESS
library(psych)

pairs.panels(data.hw, 
             method = "pearson", # Método para el cálculo de las correlaciones
             hist.col = "aquamarine",
             density = TRUE,  # Mostrando las densidades
             ellipses = FALSE #Ignorando las elipses de correlación
             )

Modelos a considerar

Modelo 1

Se considera un modelo sencillo, con la intención de obtener resultados gráficos más ilustrativos.

\(\mathbb{E}(city.mpg_i) = f(weight_i)\) donde \(city.mpg_i \thicksim normal\)

Modelo 2

\(\mathbb{E}(city.mpg_i) = f(weight_i)+ f(eng.cc_i) + f(hp_i)\) donde \(city.mpg_i \thicksim normal\)

Modelo 3

\(\mathbb{E}(hw.mpg_i) = f(weight_i)+ f(eng.cc_i) + f(hp_i)\) donde \(hw.mpg_i \thicksim normal\)

Para ajustar los modelos anteriores se utlizará un 80% como datos de entrenamiento y un 20% como datos de prueba.

train.index.city <- sample(1:nrow(data.city), 0.8 * nrow(data.city))
test.index.city <- setdiff(1:nrow(data.city), train.index.city)

train.index.hw <- sample(1:nrow(data.hw), 0.8 * nrow(data.hw))
test.index.hw <- setdiff(1:nrow(data.hw), train.index.hw)

train.city <- data.city[train.index.city,]
test.city <- data.city[test.index.city,]
  
train.hw <- data.hw[train.index.hw,]
test.hw <- data.hw[test.index.hw,]

Entrenando el modelo 1

library(mgcv)
mod1 <- gam(city.mpg ~ s(weight, k = 15), family = gaussian(link = "identity"), data = train.city)

## Por defecto la base que se utiliza es: thin plate regression splines bs = "tp".

with(train.city, plot(weight, city.mpg))

city.mpg.e <- predict.gam(mod1, data = train.city)

lines(sort(train.city$weight), predict(mod1, data = train.city)[order(train.city$weight)], col='red', type='l', lwd = 2) 

mod1$coefficients
##  (Intercept)  s(weight).1  s(weight).2  s(weight).3  s(weight).4  s(weight).5 
##  25.42073171  -4.67403012  -0.35093365   2.36434275   0.54445950  -0.07647432 
##  s(weight).6  s(weight).7  s(weight).8  s(weight).9 s(weight).10 s(weight).11 
##   0.75713648   0.01338543   0.41145302  -0.09251835  -0.51140658   0.04140351 
## s(weight).12 s(weight).13 s(weight).14 
##   0.52258899  -3.18602990 -10.95493739

Luego de entrenar el modelo se predice usando los datos de prueba.

Probando el modelo 1

city.mpg.p <- predict.gam(mod1, newdata = test.city[,-4])

# Base de datos de observados vs predecidos

header <- c("city.mpg.p", "city.mpg.o") ##Creando el nombre de las variables

datos.o.p <- data.frame(city.mpg.p, test.city[ ,4]) 

colnames(datos.o.p) <- header ##Nombrando las variables

library(tidyr)

l <- reshape(datos.o.p, 
             varying = header, 
             v.names = "mpg",
             timevar = "PredvsObs", 
             times = header, 
             direction = "long")  ##Convirtiendo a forma larga los datos 

weight.obs <- c(test.city[,1], test.city[,1]) ##Creando columna de weight

largo <- l[,-3]  ##Eliminado el indicador de la observación

datos.o.p.d <- cbind(largo,weight.obs) ##Concatenando las variables

#Construcción del gráfico de predicción vs. observado

g2 <- ggplot(datos.o.p.d, aes(x=weight.obs, y=mpg, color=PredvsObs)) +
  geom_point()

g2

En la gráfica se observa que el modelo es muy sencillo para lograr una buena predicción.

eyc <- test.city[,4] - city.mpg.p

eyc = eyc[!is.na(eyc)]

mean((eyc)^2)
## [1] 10.63697

Entrenando el modelo 2

library(mgcv)
mod2 <- gam(city.mpg ~ s(weight) + s(eng.cc) + s(hp), family = gaussian(link = "identity"), data = train.city)

city.mpg.e2 <- predict.gam(mod2, data = train.city)

mod2$coefficients
## (Intercept) s(weight).1 s(weight).2 s(weight).3 s(weight).4 s(weight).5 
##  25.4506173  -4.8882752   7.8243518  -4.1128296   3.9078352  -0.4412682 
## s(weight).6 s(weight).7 s(weight).8 s(weight).9 s(eng.cc).1 s(eng.cc).2 
##  -3.2085514  -0.2650221 -17.1523101  -7.9782309  -8.9973121  10.1427456 
## s(eng.cc).3 s(eng.cc).4 s(eng.cc).5 s(eng.cc).6 s(eng.cc).7 s(eng.cc).8 
##   5.3935129 -10.2738117   3.5453529   5.3781296  -4.1140979  12.4244031 
## s(eng.cc).9     s(hp).1     s(hp).2     s(hp).3     s(hp).4     s(hp).5 
##   6.4912470   9.0759972  11.9167114  -5.0470124   7.6893103  -4.8985615 
##     s(hp).6     s(hp).7     s(hp).8     s(hp).9 
##  -5.9026798   4.5130897 -19.4859272  -9.1702485

Probando el modelo 2

gam.check(mod2)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 11 iterations.
## The RMS GCV score gradient at convergence was 9.311271e-05 .
## The Hessian was positive definite.
## Model rank =  28 / 28 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##             k'  edf k-index p-value  
## s(weight) 9.00 6.65    1.05    0.74  
## s(eng.cc) 9.00 6.78    0.85    0.02 *
## s(hp)     9.00 6.28    1.02    0.56  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
city.mpg.p2 <- predict.gam(mod2, newdata = test.city[,-4])

eyc2 <- test.city[,4] - city.mpg.p2

eyc2 = eyc2[!is.na(eyc2)]

mean((eyc2)^2)
## [1] 4.735543

Entrenando el modelo 3

library(mgcv)
mod3 <- gam(hw.mpg ~ s(weight) + s(eng.cc) + s(hp), family = gaussian(link = "identity"), data = train.hw)

mod3$coefficients
## (Intercept) s(weight).1 s(weight).2 s(weight).3 s(weight).4 s(weight).5 
##  30.7407407  -2.5903283  -5.6966901  -1.1322530  -3.1184255   0.2187699 
## s(weight).6 s(weight).7 s(weight).8 s(weight).9 s(eng.cc).1 s(eng.cc).2 
##  -2.6225416  -0.2818628 -13.6828062  -5.4885843  -3.3857895   1.8722094 
## s(eng.cc).3 s(eng.cc).4 s(eng.cc).5 s(eng.cc).6 s(eng.cc).7 s(eng.cc).8 
##   1.3317755   2.3472869  -0.8209401  -1.8322373  -1.4682960   7.5326922 
## s(eng.cc).9     s(hp).1     s(hp).2     s(hp).3     s(hp).4     s(hp).5 
##   2.1029141   9.4474870  14.2643426  -6.1810792   9.2832579   5.2201110 
##     s(hp).6     s(hp).7     s(hp).8     s(hp).9 
##  -6.7040547   5.8467650 -26.6188147  -8.2828561

Probando el modelo 3

hw.mpg.p <- predict.gam(mod3, newdata = test.hw[,-4])

ey <- test.hw[,4] - hw.mpg.p

ey = ey[!is.na(ey)]

mse3 <- na.omit(mean((ey)^2))

mse3
## [1] 8.847609

Finalmente se compara el modelo 3 con un modelo sencillo donde se relacione la variable hw.mpg con la variable weight

library(mgcv)
mod4 <- gam(hw.mpg ~ s(weight), family = gaussian(link = "identity"), data = train.hw)

mod4$coefficients
## (Intercept) s(weight).1 s(weight).2 s(weight).3 s(weight).4 s(weight).5 
##  30.7439024  -4.5849919   3.6650385  -1.8661778  -2.5881306  -0.2635514 
## s(weight).6 s(weight).7 s(weight).8 s(weight).9 
##  -2.3517142  -0.2330649 -10.8165482 -10.4106663
hw.mpg.p1 <- predict.gam(mod4, newdata = test.hw[,-4])

ey1 <- test.hw[,4] - hw.mpg.p1

ey1 = ey1[!is.na(ey1)]

mse4 <- mean((ey1)^2)

mse4
## [1] 18.40791

Otra vez se evidencia la poca capacidad predictiva del modelo sencillo.

Modelo con distribución gamma

Finalmente, se puede comparar el resultado del modelo 2 suponiendo que \(y_i\sim gama\) con función de enlace logaritmica.

library(mgcv)
mod5 <- gam(city.mpg ~ s(weight) + s(eng.cc) + s(hp), family = Gamma(link = log), data = train.city)
mod5$coefficients
##  (Intercept)  s(weight).1  s(weight).2  s(weight).3  s(weight).4  s(weight).5 
##  3.205039027 -0.143798305  0.270903073 -0.156645569  0.139969621 -0.010281012 
##  s(weight).6  s(weight).7  s(weight).8  s(weight).9  s(eng.cc).1  s(eng.cc).2 
## -0.117948227 -0.005673164 -0.615154015 -0.277799148 -0.453115355  0.580726249 
##  s(eng.cc).3  s(eng.cc).4  s(eng.cc).5  s(eng.cc).6  s(eng.cc).7  s(eng.cc).8 
##  0.301961542 -0.589988491  0.194906265  0.312547037 -0.246432423  0.782294165 
##  s(eng.cc).9      s(hp).1      s(hp).2      s(hp).3      s(hp).4      s(hp).5 
##  0.292101961  0.244634862  0.329010441 -0.135463103  0.248215375 -0.158963234 
##      s(hp).6      s(hp).7      s(hp).8      s(hp).9 
## -0.169693205  0.132977729 -0.683801039 -0.248944231
gam.check(mod5)

## 
## Method: GCV   Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-1.260891e-10,4.524988e-10]
## (score 0.007410235 & scale 0.00655104).
## Hessian positive definite, eigenvalue range [5.554587e-05,0.0001098005].
## Model rank =  28 / 28 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##             k'  edf k-index p-value   
## s(weight) 9.00 6.60    0.99    0.45   
## s(eng.cc) 9.00 7.39    0.83    0.01 **
## s(hp)     9.00 6.28    1.01    0.52   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
city.mpg.p5 <- exp(predict.gam(mod5, newdata = test.city[,-4]))

eyc5 <- test.city[,4] - city.mpg.p5

eyc5 = eyc5[!is.na(eyc5)]

mean((eyc5)^2)
## [1] 4.651634