¿Qué son los modelos GAMLSS?

Los modelos aditivos generalizados para posición, escala y forma GAMLSS (Generalized Additive Models for Location, Scale and Shape), son modelos de regresión semi-paramétricos. Paramétricos en cuanto a que requieren asumir que la variable respuesta sigue una determinada distribución paramétrica (normal, beta, gamma…) y semi porque los parámetros de esta distribución pueden ser modelados, cada uno de forma independiente, siguiendo funciones no paramétricas (lineales, aditivas o no lineales). Esta versatilidad hace de los GAMLSS una herramienta adecuada para modelar variables que siguen todo un abanico de distribuciones (no normales, asimétricas, con varianza no constante…).

Aunque en este documento se intentan resumir los principales aspectos prácticos, son muchos los aspectos teóricos (distribuciones, métodos de ajuste, penalizaciones…) que engloban de este tipo de modelos. Para profundizar en ellos, la mejor opción es el libro escrito por los propios autores de los modelos GAMLSS.

Stasinopoulos, Dm & Rigby, Robert & Heller, Gillian & Voudouris, Vlasios & De Bastiani, Fernanda. (2017). Flexible regression and smoothing: Using GAMLSS in R. 10.1201/b21973.

¿Cómo funcionan?

Los modelos GAMLSS asumen que la variable respuesta tiene una función de densidad definida por hasta 4 parámetros \((\mu, \sigma, \nu, \tau)\) que determinan su posición (p.ej. media), escala (p.ej. desviación estándar) y forma (p. ej. skewness y kurtosis), y que cada uno de ellos puede variar independientemente de los otros en función de los predictores. Estos modelos aprenden por lo tanto hasta 4 funciones, donde cada una establece la relación entre las variables predictoras y uno de los parámetros.

Son muchas las distribución que pueden emplearse para la variable respuesta \(y\), \(f(y_i|\mu_i, \sigma_i, \nu_i, \tau_i)\), la única condición es que el logaritmo de la función de densidad (\(\log f(y_i|\mu_i, \sigma_i, \nu_i, \tau_i)\)) y su primera derivada respecto a cada uno de los parámetros puedan calcularse.

Ventajas de los modelos GAMLSS

Para entender mejor las ventajas de este tipo de modelos, conviene repasar brevemente la evolución que han seguido los modelos estadísticos desde el modelo lineal (LM), modelo lineal generalizado (GLM), modelo aditivo generalizado (GAM) y finalmente los modelos GAMLSS.

Modelos de regresión lineal (LM)

El modelo de regresión lineal (Legendre, Gauss, Galton y Pearson) considera que, dado un conjunto de observaciones, la media \(\mu\) de la variable respuesta \(Y\) se relaciona de forma lineal con la o las variables regresoras \(X\).

\[\mu = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}\]

o en notación matricial (incorporando \(\beta_0\) en el vector \(\mathbf{\beta}\):

\[\mu = \mathbf{X}^T \mathbf{\beta}\]

Ajustar el modelo consiste en estimar los valores de los coeficientes de regresión \(\hat{\beta}\) y la varianza \(\hat{\sigma}^2\) que maximizan la verosimilitud (likelihood) de los datos, es decir, los que dan lugar al modelo que con mayor probabilidad puede haber generado los datos observados. El método empleado con más frecuencia es el ajuste por mínimos cuadrados (OLS):

\[\hat{\beta} = \underset{\beta_0,\beta}{\arg\min} (Y - \mathbf{X}^T \mathbf{\beta})^2\]

\[\hat{\mu} = \mathbf{X}^T \mathbf{\hat{\beta}}\]

\[\hat{\sigma}^2 = \frac{\sum^n_{i=1} \hat{\epsilon}_i^2}{n-p} = \frac{\sum^n_{i=1} (y_i - \hat{\mu})^2}{n-p}\]

donde \(\hat{\epsilon}_i\) es la diferencia entre cada valor observado y la media estimada, \(n\) es el número de observaciones y \(p\) el número de predicciones.

Es importante resaltar una característica de estos modelos. La estimación de la media \(\hat{\mu}\) depende directamente del valor de los predictores \(\mathbf{X}\), sin embargo, la estimación de la varianza \(\hat{\sigma}^2\) no depende directamente de los predictores, sino de \(\hat{\mu}\). Para que sean válidas, deben cumplirse varias condiciones, dos de las cuales son:

• La variable respuesta \(Y\) se tiene que distribuir de forma normal.

• La varianza \(\sigma^2\) debe de ser constante en todo el rango de la variable respuesta.
  
  

Modelos lineales generalizados (GLM)

El requisito de que la variable respuesta se distribuya de forma normal hace que el modelo de regresión lineal no sea aplicable a muchos problemas reales. Los modelos GLM (John Nelder y Robert Wedderburn) permiten superar esta limitación haciendo una abstracción a un modelo más genérico, en el que la variable respuesta puede seguir cualquier distribución de la familia exponencial (Gaussian, Poisson, binomial…). Para conseguirlo, la media de la variable respuesta ya no está relacionada directamente con los predictores, sino que se hace a través de una función link \(g(.)\). Por ejemplo, el modelo lineal generalizado es equivale al modelo lineal por mínimos cuadrados cuando la función link es la identidad, o a la regresión logística cuando es binomial (logit).

\[g(\mu) = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}\]

\[g(\mu) = \mathbf{X}^T \mathbf{\beta}\]

Cuando se emplea una función identidad, la máxima verosimilitud equivale a minimizar la suma de errores cuadrados, lo que tiene una solución analítica. Para el resto de familias o cuando se emplea la regularización L1, no existe una solución analítica, por lo que se requiere un método de optimización iterativo como el iteratively reweighted least squares (IRLSM), L-BFGS, método de Newton o descenso de gradiente.

Modelos de regresión lineal generalizados aditivos (GAM)

A pesar de las posibilidades que ofrecen los GLM, siguen teniendo como limitación que la relación entre los predictores y la media de la variable respuesta debe de ser lineal y constante. En 1990 Hastie y Tibshirani introdujeron los modelos GAM, una extensión de los modelos GLM que permite incorporar relaciones no lineales. En los modelos GAM, la relación de cada predictor \(x_i\) con la media de la variable respuesta \(g(\mu)\) no es directa, sino que se hace a través de una función \(f(x_i)\).

\[\eta = g(\mu) = \beta_0 + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\]

La función \(f_i(x_i)\) puede ser cualquier función que recibe el valor del predictor \(x_i\) y devuelve otro valor, pudiendo ser esta transformación lineal o no lineal. En la práctica, las funciones más empleadas son funciones no lineales tipo smooth: cubic regression splines, thin plate regression splines o Penalized splines.

Modelos aditivos generalizados para posición, escala y forma (GAMLSS)

Los modelos lineales generalizados (GLM) y los modelos generalizados aditivos (GAM) asumen que la variable respuesta \(Y\) sigue una distribución de la familia exponencial, cuya media \(\mu\) puede ser modelada en función de otras variables (predictores) y cuya varianza \(\sigma\) se calcula mediante una constante de dispersión y una función \(\nu(\mu)\). Esto último significa que, la varianza, skewness y kurtosis, no se modelan directamente en función de las variables predictoras sino de forma indirecta a través de su relación con la media. ¿Qué impacto puede tener esto en la práctica? Si bien las estimación de la media obtenidas por estos modelos son buenas, no lo son tanto las incertidumbres asociadas y, en consecuencia, los intervalos de confianza y de predicción que se puedan calcular.

Los modelos GAMLSS, introducidos por Rigby y Stasinopoulos en 2005, permiten, además de incorporar distribuciones que no son de la familia exponencial, modelar explícitamente cada uno de sus parámetros en función de las variables predictoras empleando funciones lineales y no lineales. Los términos empleados dentro del marco de los GAMLSS para referirse a los parámetros de localización, escala y forma son \((\mu, \sigma, \nu, \tau)\) .

\[\mathbf{Y} \sim D(\mu, \sigma, \nu, \tau)\]

\[\mathbf{Y} = \mathbf{X}^T \mathbf{\beta}\]

\[\eta_1 = g_1(\mu) = \mathbf{X}^T \mathbf{\beta} + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\]

\[\eta_2 = g_2(\sigma) = \mathbf{X}^T \mathbf{\beta} + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\]

\[\eta_3 = g_3(\nu) = \mathbf{X}^T \mathbf{\beta} + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\]

\[\eta_4 = g_4(\tau) = \mathbf{X}^T \mathbf{\beta} + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\] donde \(\mathbf{Y} \sim D(\mu, \sigma, \nu, \tau)\) es la distribución de la variable respuesta (pueden ser menos parámetros), \(\mathbf{X}\) contiene los términos lineales del modelo, \(\mathbf{\beta}\) son los coeficientes lineales y \(f_i(x_i)\) son funciones de suavizado no lineales (smooth) de cada predictor.

Los GAMLSS son una forma de superar las limitaciones de los modelos GLM y GAM. Como resultado del modelo, se consigue caracterizar la distribución completa, permitiendo generar intervalos probabilísticos y predicción de cuantiles.

Ajuste de los modelos GAMLSS

Los modelos GAMLSS se ajustan mediante una adaptación del algoritmo backfitting, un algoritmo típicamente empleado para el ajuste de modelos aditivos. Por ejemplo, para el modelo aditivo

\[g(\mu) = \beta_0 + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\]

en términos generales, el algoritmo de ajuste por backfitting funciona de la siguiente forma:

1. Se inicia con un valor para todos los términos \(f_i\) del modelo, por ejemplo poniéndolos todos a cero.

2. Se estima el valor del primer término \(f_1\) empleando los datos de entrenamiento, con el objetivo de predecir lo mejor posible la variable respuesta \(y\).

3. Se estima el valor del segundo término \(f_2\) empleando residuos \(y - f_1(x_1)\).

4. Se repite el paso 3, ajustando cada término con los residuos del anterior.

5. Tras ajustar todos los términos, el valor del primer término se descarta y se reajusta empleando el residuo de todos los demás términos.

6. Repetir todos los pasos del 3 al 5 hasta que se alcance un criterio de parada (que los valores apenas cambien o que se alcance un número máximo de iteraciones).

Selección de modelos GAMLSS

Para identificar el mejor modelo GAMLSS es necesario comparar diferentes modelos candidatos, cada uno con una combinación de distribución para la variable respuesta, función link para cada uno de los parámetros, predictores e hiperparámetros. Existen varias estrategias para la selección del mejor modelo de entre los comparados:

• Generalized Akaike information criterion (GAIC): esta métrica emplea el log likelihood del modelo multiplicado por \(-2\) y añade una penalización por cada parámetro que incluye el modelo. Más detalles sobre el GAIC en el apartado métricas de ajuste.

• Generalized cross-validation y cross-validation

Acorde a la comunidad de estadísticos, el segundo método es más recomendable, aunque supone mayor coste computacional.

gamlss

La implementación en R de los modelos GAMLSS está disponible en el paquete gamlss y todas sus extensiones:

• gamlss: paquete central con las principales funciones para ajustar modelos GAMLSS.

• gamlss.dist: extensión con distribuciones adicionales.

• gamlss.cens: extensión para variables respuesta censuradas.

• gamlss.mx: extensión para distribuciones mixtas (combinación de distribuciones).

• gamlss.nl: extensión para ajustar modelos no lineales.

• gamlss.tr: extensión para distribuciones truncadas.

• gamlss.add: extensión que permite incluir nuevas transformaciones no lineales a los predictores de un modelo aditivo (árboles, redes neuronales y splines multidimensionales).

La principal función del paquete es gamlss(), con la que se ajustan y generan los modelos. A lo largo de este documento se muestran, mediante ejemplos prácticos, como emplear estos paquetes para crear modelos GAMLSS.

Datos

El set de datos rent del paquete gamlss.data contiene información sobre la precio de 1969 viviendas situadas en Munich en el año 1993. Además del precio, incluye 9 variables adicionales:

• R: precio del alquiler

• Fl: metros cuadrados de la vivienda

• A: año de construcción

• Sp: si la calidad del barrio donde está situada la vivienda es superior la media (1) o no (0).

• Sm: si la calidad del barrio donde está situada la vivienda es inferior la media (1) o no (0).

• B: si tiene cuarto de baño (1) o no (0).

• H: si tiene calefacción central (1) o no (0).

• L: si el equipamiento de la cocina está por encima de la media (1) o no (0).

• loc: combinación de Sp y Sm indicando si la calidad del barrio donde está situada la vivienda es inferior (1), igual (2) o superior (3) a la media.

El objetivo es obtener un modelo capaz de predecir el precio del alquiler. Siguiendo el ejemplo de Stasinopoulos, Dm & Rigby, Robert & Heller, Gillian & Voudouris, Vlasios & De Bastiani, Fernanda , se emplean como predictores únicamente Fl, A, H y loc.

library(gamlss)
library(tidyverse)
library(ggpubr)
library(skimr)

data("rent")
datos <- rent
datos <- datos %>% select(R, Fl, A, H, loc)
# Se renombran las variables para que sean más explicativas
colnames(datos) <- c("precio", "metros", "anyo", "calefaccion", "situacion")
skim(datos)

Variable type: factor

Variable type: numeric

Los siguientes gráficos muestran la distribución del precio de las viviendas y su relación con cada uno de los predictores.

p1 <- ggplot(data = datos, aes(x = metros, y = precio)) + 
      geom_point(alpha = 0.4) +
      labs(ttile = "Precio vs metros") +
      theme_bw()

p2 <- ggplot(data = datos, aes(x = anyo, y = precio)) + 
      geom_point(alpha = 0.4) +
      labs(ttile = "Precio vs año") +
      theme_bw()

p3 <- ggplot(data = datos, aes(x = calefaccion, y = precio)) + 
      geom_violin() +
      geom_boxplot(width = 0.1, outlier.shape = NA) +
      labs(ttile = "Precio vs calefacción") +
      theme_bw()

p4 <- ggplot(data = datos, aes(x = situacion, y = precio)) + 
      geom_violin() +
      geom_boxplot(width = 0.1, outlier.shape = NA) +
      labs(ttile = "Precio vs situación") +
      theme_bw()

ggpubr::ggarrange(
  plotlist = list(p1, p2, p3, p4)
) %>%
ggpubr::annotate_figure(
  top = text_grob("Relación entre el precio y el resto de variables",
                  color = "Black",
                  face  = "bold",
                  size  = 14,
                  x     = 0.3)
  )

ggplot(data = datos, aes(x = precio)) +
  geom_density(alpha = 0.5, fill = "gray50") +
  geom_rug(alpha = 0.2) +
  labs(title = "Distribución del precio de los pisos") +
  theme_bw()

Con estos gráficos se evidencia que:

• La relación entre el precio y los metros cuadrados es lineal y positiva, con una varianza creciente a medida que aumentan los metros cuadrados. Esto significa que no se puede asumir la condición de varianza constante (homocedasticidad).

• La distribución del precio es asimétrica con una cola positiva, la variable respuesta no se distribuye de forma normal.

• El precio medio de pisos con calefacción es superior a los que no tienen calefacción.

• El precio de los pisos aumenta progresivamente si está en un barrio por debajo de la media, en la media o por encima.

Dadas estas características, se necesita un modelo que:

• Sea capaz de aprender relaciones complejas entre varios predictores, incluyendo relaciones no lineales.

• Sea capaz de modelar explícitamente la varianza en función de los predictores, ya que esta no es constante.

• Sea capaz de aprender distribuciones asimétricas con una marcada cola positiva.
  
  

Regresión Lineal (LM)

# El resultado es igual al obtenido con lm()
modelo_lm <- gamlss(
              formula = precio ~ metros + anyo + calefaccion + situacion,
              family  = NO,
              data    = datos,
              trace   = FALSE
             )

summary(modelo_lm)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = precio ~ metros + anyo + calefaccion +  
##     situacion, family = NO, data = datos, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2775.0388   470.1352  -5.903 4.20e-09 ***
## metros           8.8394     0.3370  26.228  < 2e-16 ***
## anyo             1.4808     0.2385   6.208 6.55e-10 ***
## calefaccion1  -204.7596    18.9858 -10.785  < 2e-16 ***
## situacion2     134.0523    25.1430   5.332 1.09e-07 ***
## situacion3     209.5815    27.1286   7.725 1.76e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.73165    0.01594   359.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1969 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  1962 
##                       at cycle:  2 
##  
## Global Deviance:     28159 
##             AIC:     28173 
##             SBC:     28212.1 
## ******************************************************************

La función link para estimar la varianza es \(\log\), ya que la varianza no puede ser negativa. Para obtener el valor estimado de la varianza en la misma escala que la variable respuesta, hay que aplicar el exponente.

paste("El valor estimado de la varianza es:", exp(5.73165))

## [1] "El valor estimado de la varianza es: 308.477837124068"

Con la función plot() se muestra el diagnostico de los residuos del modelo.

plot(modelo_lm)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  4.959554e-13 
##                        variance   =  1.000508 
##                coef. of skewness  =  0.7470097 
##                coef. of kurtosis  =  4.844416 
## Filliben correlation coefficient  =  0.9859819 
## ******************************************************************

Los gráficos worm son otra forma de evaluar, visualmente, la calidad de un modelo a través de sus residuos. Esta representación es similar a la de un gráfico Q-Q pero sustrayendo la línea teórica, por lo que, idealmente, los valores resultantes deberían de ser cero (línea horizontal discontinua). Las dos curvas elípticas discontinuas muestran el intervalo de confianza del 95%. Si el modelo es correcto, sólo entorno al 5% de las observaciones deberían quedar fuera. Por defecto, la función wp() muestra también un ajuste cúbico (curva continua roja) que ayuda a identificar la tendencia de los residuos.

# Worm plot de los residuos
wp(modelo_lm, ylim.all = 1)

El estudio de los residuos muestra que:

• Los errores no se distribuyen de forma aleatoria en torno a cero, sino que aumentan a medida que aumenta el valor predicho. Esto es indicativo y consecuencia de la heterogeneidad de la varianza, en concreto, la varianza aumenta con la media. (Gráfico superior-izquierda)

• La condición de normalidad no se cumple, el gráfico Q-Q y el gráfico worm muestran una clara desviación entre los cuantiles observados y los teóricos.

La violación de la condición de distribución normal conlleva que el modelo LM no sea adecuado para captar la relación entre el precio y los predictores.

Lineal Generalizado (GLM)

La distribución de la variable precio tiene dos características: no puede ser negativa y tiene una cola positiva (asimetría). Esto hace de la distribución gamma una buena candidata. Se puede encontrar un listado con todas las distribuciones disponibles en gamlss en este documento.

La distribución gamma tiene dos parámetros que tienen que ser estimados: media (\(\mu\)) y escala (\(\sigma\)). Por defecto, el paquete gamlss emplea para ambos la función \(log\) como link, ya que los dos parámetros tienen un dominio (0, \(\infty\)).

\[\eta_1 = g_1(\mu) = \log(\mu)\]

\[\eta_2 = g_1(\sigma) = \log(\sigma)\]

Es muy importante tener en cuenta este aspecto a la hora de utilizar gamlss, ya que, dependiendo de qué función link que se emplee, una misma distribución puede dar lugar a modelos distintos. Por defecto, este paquete emplea funciones link que reflejan el dominio del parámetro: “identity” para (\(-\infty\), \(\infty\)), “log” para (\(0\), \(\infty\)), “logit” para (\(0\), \(1\)), etc.

# Mostrar las funciones link por defecto de la distribución gamma.
GA()

## 
## GAMLSS Family: GA Gamma 
## Link function for mu   : log 
## Link function for sigma: log

# Mostrar las funciones link disponibles para la distribución gamma.
show.link(family = GA)

## $mu
## c("inverse", "log", "identity", "own")
## 
## $sigma
## c("inverse", "log", "identity", "own")

modelo_glm <- gamlss(
                formula = precio ~ metros + anyo + calefaccion + situacion,
                family  = GA,
                data    = datos,
                trace   = FALSE
              )

summary(modelo_glm)

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  gamlss(formula = precio ~ metros + anyo + calefaccion +  
##     situacion, family = GA, data = datos, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.8649770  0.5681265   5.043 5.01e-07 ***
## metros        0.0106232  0.0004128  25.735  < 2e-16 ***
## anyo          0.0015100  0.0002886   5.232 1.86e-07 ***
## calefaccion1 -0.3000745  0.0231141 -12.982  < 2e-16 ***
## situacion2    0.1907641  0.0305197   6.251 5.00e-10 ***
## situacion3    0.2640828  0.0329216   8.022 1.78e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.98220    0.01558  -63.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1969 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  1962 
##                       at cycle:  2 
##  
## Global Deviance:     27764.59 
##             AIC:     27778.59 
##             SBC:     27817.69 
## ******************************************************************

El modelo resultante es:

\[\text{precio} \sim GA(\hat{\mu}, \hat{\sigma})\]

donde:

\(log(\hat{\mu}) = 2.865 + 0.0106\text{metros} + 0.00151\text{año} − \\0.3 \ (\text{si calefaccion=1}) + 0.191 \ (\text{si localicacion=2}) + 0.264 \ (\text{si localicacion=3})\)

\(log(\hat{\sigma}) = -0.9822\)

plot(modelo_glm)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  0.0004795675 
##                        variance   =  1.000657 
##                coef. of skewness  =  -0.1079453 
##                coef. of kurtosis  =  3.255464 
## Filliben correlation coefficient  =  0.9990857 
## ******************************************************************

# Worm plot de los residuos
wp(modelo_glm, ylim.all = 0.5)

El análisis de los residuos muestra un diagnóstico más favorable que en el caso del modelo lineal LM. El gráfico de los residuos vs el valor real (esquina superior izquierda) tiene menor heterogeneidad, y la curvatura en las colas del gráfico Q-Q (esquina inferior derecha) se ha reducido notablemente. Esto apunta a que el modelo GLM con una distribución gamma es más adecuado para modelar el precio de las viviendas. ¿Cómo corroborarlo? Una forma de hacerlo es comparando la métrica GAIC de ambos modelos.

# Por defecto, GAIC emplea una penalización k = 2.5
GAIC(modelo_lm, modelo_glm)

El modelo GLM consigue reducir el valor GAIC. Teniendo en cuenta esto, junto con el diagnóstico de los residuos, puede considerarse mejor candidato que el modelo LM.

Generalizado Aditivo (GAM)

En los modelos GAM se pueden aplicar funciones (lineales o no lineales) a cada uno de los predictores. El paquete gamlss implementa varias funciones smooth que suelen dar buenos resultados. En este caso se emplea P-splines (pb()).

# Se emplea P-splines para los predictores continuos
# Distribución gamma para la variable respuesta
modelo_gam <- gamlss(
                formula = precio ~ pb(metros) + pb(anyo) + calefaccion + situacion,
                family  = GA,
                data    = datos,
                trace   = FALSE
             )

summary(modelo_gam)

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  gamlss(formula = precio ~ pb(metros) + pb(anyo) + calefaccion +  
##     situacion, family = GA, data = datos, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.0851197  0.5691903   5.420 6.69e-08 ***
## pb(metros)    0.0103084  0.0004031  25.574  < 2e-16 ***
## pb(anyo)      0.0014062  0.0002892   4.862 1.26e-06 ***
## calefaccion1 -0.3008111  0.0225874 -13.318  < 2e-16 ***
## situacion2    0.1886692  0.0299282   6.304 3.57e-10 ***
## situacion3    0.2719856  0.0322843   8.425  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.00196    0.01559  -64.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1969 
## Degrees of Freedom for the fit:  11.21547
##       Residual Deg. of Freedom:  1957.785 
##                       at cycle:  3 
##  
## Global Deviance:     27683.22 
##             AIC:     27705.65 
##             SBC:     27768.29 
## ******************************************************************

Debido a la incorporación de las funciones no lineales (smoothers), hay que ser cauto a la hora de interpretar los coeficientes y los errores mostrados en el summary. Los coeficientes y errores de los predictores con funciones smooth (metros y anyo) contemplan únicamente la parte lineal, ignorando la no lineal. En el caso de los predictores lineales (calefacción y situacion) sus errores se estiman asumiendo que los términos con smoother son fijos, y no contemplan la incertidumbre introducida al ajustar cada una de las funciones smooth.

Para conocer la contribución total (lineal y no lineal) de los predictores transformados por funciones smooth se puede emplear la función drop1(). Esta función calcula el impacto que tiene en el modelo (en términos de grados de libertad totales, AIC y significancia estadística) el eliminar cada uno de los predictores de forma secuencial.

# Esta función puede tardar si hay muchos predictores no lineales o muchos datos.
drop1(modelo_gam, parallel = "multicore", ncpus = 4)

Otra consecuencia de añadir términos no lineales en un modelo es que su descripción matemática no es fácilmente interpretable, en su lugar, es más ilustrativo representar la relación de cada predictor con la variable respuesta. Esto puede hacerse con la función term.plot().

term.plot(modelo_gam, pages = 1, ask = FALSE, rug = TRUE)

Los gráficos muestran la relación entre \(\eta = log(\mu)\) (el logaritmo de la media del precio de las viviendas) con cada uno de los predictores:

• El valor aumenta linealmente a medida que aumentan los metros de la vivienda.

• La relación con el anyo de construcción no es lineal: se mantiene constante hasta aproximadamente 1960 y luego aumenta.

• La interpretación de los dos otros predictores (calefaccion y situación) es la misma que en los modelos anteriores, ya que se han incorporado al modelo sin sufrir ninguna transformación.

plot(modelo_gam)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  0.0005312226 
##                        variance   =  1.00067 
##                coef. of skewness  =  -0.1235641 
##                coef. of kurtosis  =  3.322109 
## Filliben correlation coefficient  =  0.9989678 
## ******************************************************************

# Worm plot de los residuos
wp(modelo_gam, ylim.all = 0.5)

Se compara el modelo GAM con los dos anteriores:

GAIC(modelo_lm, modelo_glm, modelo_gam)

El modelo GAM ha conseguido reducir todavía más el valor AIC y la distribución de sus residuos ha mejorado ligeramente.

Generalizado Aditivo Localización, Escala y Forma (GAMLSS)

Existen 3 formas distintas de parametrizar la distribución gamma, una de ellas es mediante la media \(\mu\) y un parámetro de escala \(\theta\) (llamado \(\sigma\) en el paquete gamlss). Hasta el momento, se ha modelado únicamente \(\mu\) en función de los predictores, asumiendo que la escala es constante y puede ser estimada a partir de la media.

Una forma de mejorar el modelo es modelando también su parámetro de escala en función de los predictores. Esto resulta particularmente útil cuando los datos muestran heterocedasticidad (su varianza no es constante). Se repite el ajuste, pero esta vez modelando sus dos parámetros en función de los predictores.

\[Y \sim GA(\mu, \sigma)\]

\[\eta_1 = g_1(\mu) = \beta_0 + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\]

\[\eta_2 = g_2(\sigma) = \beta_0 + f_1(x_{1}) + f_2(x_{2}) + ... + f_p(x_{p})\]

# Modelo GAMLSS para distribución gamma con parámetros de media y escala
# Se emplea P-splines para los predictores continuos
modelo_gamlss <- gamlss(
                  formula = precio ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                  sigma.formula = ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                  family  = GA,
                  data    = datos,
                  trace   = FALSE
                )

summary(modelo_gamlss)

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  gamlss(formula = precio ~ pb(metros) + pb(anyo) + calefaccion +  
##     situacion, sigma.formula = ~pb(metros) + pb(anyo) +  
##     calefaccion + situacion, family = GA, data = datos,      trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.8844571  0.5843467   4.936 8.64e-07 ***
## pb(metros)    0.0105402  0.0003620  29.116  < 2e-16 ***
## pb(anyo)      0.0014977  0.0002966   5.050 4.83e-07 ***
## calefaccion1 -0.2918923  0.0242854 -12.019  < 2e-16 ***
## situacion2    0.1938688  0.0323332   5.996 2.40e-09 ***
## situacion3    0.2734326  0.0339197   8.061 1.31e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.9225569  0.8645448   6.850 9.84e-12 ***
## pb(metros)    0.0019227  0.0007734   2.486  0.01300 *  
## pb(anyo)     -0.0035795  0.0004372  -8.187 4.79e-16 ***
## calefaccion1  0.0659289  0.0415210   1.588  0.11248    
## situacion2   -0.1166325  0.0565969  -2.061  0.03946 *  
## situacion3   -0.1702147  0.0608990  -2.795  0.00524 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1969 
## Degrees of Freedom for the fit:  22.25035
##       Residual Deg. of Freedom:  1946.75 
##                       at cycle:  4 
##  
## Global Deviance:     27570.28 
##             AIC:     27614.78 
##             SBC:     27739.06 
## ******************************************************************

Al igual que se hizo anteriormente con \(\mu\), se puede evaluar la influencia de cada predictor sobre el parámetro \(\sigma\).

term.plot(modelo_gamlss, parameter = "sigma", pages = 1, ask = FALSE, rug = TRUE)

drop1(modelo_gamlss, parameter = "sigma", parallel = "multicore", ncpus = 4)

Viendo el gráfico y el impacto que tiene eliminar el predictor calefaccion del modelo, puede decirse que, en presencia de los otros, no aporta información para modelar el parámetro \(\sigma\).

modelo_gamlss <- gamlss(
                  formula = precio ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                  sigma.formula = ~ pb(metros)+pb(anyo)+situacion,
                  family  = GA,
                  data    = datos,
                  trace   = FALSE
                )

Véanse ahora los residuos.

# Worm plot de los residuos
wp(modelo_gamlss, ylim.all = 0.5)

Un estudio visual del gráfico worm indica que, con este modelo, el número de observaciones cuyos residuos están fuera del intervalo de confianza del 95% es el que cabría esperar. Aun así, la clara tendencia en U, refleja que el modelo no es lo suficientemente flexible para capturar la asimetría (cola) de la distribución.

Por último, se comparan los 3 modelos ajustados:

GAIC(
  modelo_lm,
  modelo_glm,
  modelo_gam,
  modelo_gamlss
)

Acorde al criterio GAIC, el modelo GAMLSS es el que mejor captura la relación entre el precio de la vivienda y los predictores empleados.

La distribución gamma empleada tiene únicamente dos parámetros, uno de posición \(\mu\) y otro de escala \(\sigma\). La limitación de las distribuciones con solo dos parámetros es que asumen que la asimetría (skewness y kurtosis) es constante para un valor fijo de \(\mu\) y \(\sigma\). Cuando se dispone de una cantidad elevada de datos se puede tratar de emplear distribuciones en las que sus parámetros de asimetría, llamados en este paquete \(\eta\) y \(\nu\), pueden ser modelados también en función de los predictores. Es aquí donde reside el potencial de estos modelos.

Véase el resultado de emplear una distribución con 3 parámetros Box-Cox Cole-Green o Box-Cox normal, que resulta útil para modelar datos con asimetría positiva.

BCCGo()

## 
## GAMLSS Family: BCCGo Box-Cox-Cole-Green-orig. 
## Link function for mu   : log 
## Link function for sigma: log 
## Link function for nu   : identity

# Modelo GAMLSS para distribución Box-Cox Cole-Green con parámetros de media,
# escala y forma. Se emplea P-splines para los predictores continuos
modelo_gamlss_2 <- gamlss(
                    formula = precio ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                    sigma.formula = ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                    nu.formula    = ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                    family  = BCCGo,
                    data    = datos,
                    trace   = FALSE
                  )

summary(modelo_gamlss_2)

## ******************************************************************
## Family:  c("BCCGo", "Box-Cox-Cole-Green-orig.") 
## 
## Call:  gamlss(formula = precio ~ pb(metros) + pb(anyo) + calefaccion +  
##     situacion, sigma.formula = ~pb(metros) + pb(anyo) +  
##     calefaccion + situacion, nu.formula = ~pb(metros) +  
##     pb(anyo) + calefaccion + situacion, family = BCCGo,  
##     data = datos, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.0285797  0.5869375   3.456  0.00056 ***
## pb(metros)    0.0104079  0.0003755  27.721  < 2e-16 ***
## pb(anyo)      0.0019277  0.0002973   6.485 1.12e-10 ***
## calefaccion1 -0.3213906  0.0268845 -11.955  < 2e-16 ***
## situacion2    0.1853274  0.0347087   5.340 1.04e-07 ***
## situacion3    0.2742129  0.0364029   7.533 7.57e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.6534017  0.8661933   7.681 2.48e-14 ***
## pb(metros)    0.0018721  0.0008059   2.323   0.0203 *  
## pb(anyo)     -0.0039673  0.0004385  -9.047  < 2e-16 ***
## calefaccion1  0.0819273  0.0430147   1.905   0.0570 .  
## situacion2   -0.0851622  0.0627391  -1.357   0.1748    
## situacion3   -0.1410790  0.0675848  -2.087   0.0370 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -3.160132   3.203110  -0.987    0.324  
## pb(metros)   -0.000983   0.002295  -0.428    0.669  
## pb(anyo)      0.001989   0.001619   1.229    0.219  
## calefaccion1 -0.286663   0.122130  -2.347    0.019 *
## situacion2   -0.171143   0.183091  -0.935    0.350  
## situacion3   -0.084593   0.201476  -0.420    0.675  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1969 
## Degrees of Freedom for the fit:  28.41391
##       Residual Deg. of Freedom:  1940.586 
##                       at cycle:  7 
##  
## Global Deviance:     27551.32 
##             AIC:     27608.15 
##             SBC:     27766.85 
## ******************************************************************

term.plot(modelo_gamlss_2, parameter = "nu", pages = 1, ask = FALSE, rug = TRUE)

drop1(modelo_gamlss_2, parameter = "nu", parallel = "multicore", ncpus = 4)

El estudio de los coeficientes apunta a que, únicamente el predictor calefacción, contribuye a modelar el parámetro \(\nu\).

modelo_gamlss_2 <- gamlss(
                    formula = precio ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                    sigma.formula = ~ pb(metros)+pb(anyo)+calefaccion+situacion,
                    nu.formula    = ~ calefaccion,
                    family  = BCCGo,
                    data    = datos,
                    trace   = FALSE
                  )

wp(modelo_gamlss_2)

La distribución de los residuos ha mejorado respecto a la del modelo gamma. Por último, se comprueba si la métrica AIC también es mejor.

GAIC(
  modelo_lm,
  modelo_glm,
  modelo_gam,
  modelo_gamlss,
  modelo_gamlss_2
)

A lo largo de este ejemplo se han mostrado las ventajas derivadas de la flexibilidad que ofrecen los modelos GAMLSS. En los siguientes apartados, se profundiza en más aspectos teóricos y prácticos de este tipo de modelos.

Comparación y selección

La función fitDist() ajusta toda las distribuciones paramétricas disponibles de una determinada familia, y las compara acorde al GAIC (generalized Akaike information criterion). La familia de distribuciones se especifica con el argumento type y puede ser: "realAll", "realline", "realplus", "real0to1", "counts" y "binom".

• realAll: distribuciones de la familia realline + realplus.

• realline: distribuciones continuas en el dominio \((-\infty, \infty)\): NO, GU, RG ,LO, NET, TF, TF2, PE,PE2, SN1, SN2, exGAUS, SHASH, SHASHo, SHASHo2, EGB2, JSU, JSUo, SEP1, SEP2, SEP3, SEP4, ST1, ST2, ST3, ST4, ST5, SST, GT.

• realplus: distribuciones continuas en el dominio \((0, \infty]\): EXP, GA,IG,LOGNO, LOGNO2,WEI, WEI2, WEI3, IGAMMA,PARETO2, PARETO2o, GP, BCCG, BCCGo, exGAUS, GG, GIG, LNO,BCTo, BCT, BCPEo, BCPE, GB2.

• real0to1: distribuciones continuas en el dominio \([0, 1]\): BE, BEo, BEINF0, BEINF1, BEOI, BEZI, BEINF, GB1.

• counts: distribuciones para cuentas: PO, GEOM, GEOMo,LG, YULE, ZIPF, WARING, GPO, DPO, BNB, NBF,NBI, NBII, PIG, ZIP,ZIP2, ZAP, ZALG, DEL, ZAZIPF, SI, SICHEL,ZANBI, ZAPIG, ZINBI, ZIPIG, ZINBF, ZABNB, ZASICHEL, ZINBF, ZIBNB, ZISICHEL.

• binom: distribuciones para datos binomiales: BI, BB, DB, ZIBI, ZIBB, ZABI, ZABB.
  
  

Otra alternativa disponible es la función fitDistPred(), que ajusta las distribuciones igual que fitDist() pero, en lugar de compararlas por el GAIC, emplea un conjunto de test para calcular la bondad de ajuste (global deviance).

distribuciones <- fitDist(
                    datos$precio,
                    k = 2, # esta penalización equivale al AIC
                    type = "realplus",
                    trace = FALSE,
                    try.gamlss = TRUE
                  )

## Error in solve.default(oout$hessian) : 
##   Lapack routine dgesv: system is exactly singular: U[3,3] = 0
## Error in solve.default(oout$hessian) : 
##   Lapack routine dgesv: system is exactly singular: U[4,4] = 0

distribuciones$fits %>%
  enframe(name = "distribucion", value = "GAIC") %>%
  arrange(GAIC)

El objeto devuelto por fitDist() almacena la mejor de entre todas las distribuciones probadas.

summary(distribuciones)

## *******************************************************************
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## 
## Call:  gamlssML(formula = y, family = DIST[i]) 
## 
## Fitting method: "nlminb" 
## 
## 
## Coefficient(s):
##              Estimate  Std. Error   t value   Pr(>|t|)    
## eta.mu    749.2743664   8.5740894  87.38822 < 2.22e-16 ***
## eta.sigma  -0.7520409   0.0163055 -46.12189 < 2.22e-16 ***
## eta.nu      0.2530936   0.0379940   6.66141 2.7122e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   1966 
## Global Deviance:     28607.7 
##             AIC:     28613.7 
##             SBC:     28630.4

Simulación

Todas las distribuciones de los paquetes gamlss.* disponen de funciones d, p, q y r para calcular probabilidad, densidad, cuantiles y generar valores aleatorios.

BCCG()

## 
## GAMLSS Family: BCCG Box-Cox-Cole-Green 
## Link function for mu   : identity 
## Link function for sigma: log 
## Link function for nu   : identity

# La función link de BCCG() es el log, por lo que el valor devuelto por el
# ajuste está en logaritmo
rBCCG(n = 1, mu = 749.2743664, sigma = exp(-0.7520409), nu = 0.2530936)

## [1] 1503.38

Cuando se selecciona una distribución, conviene simular nuevos valores para ver si estos se parecen a los datos observados. Empleando la distribución Box-Cox-Cole-Green (BCCG) con los parámetros estimados por fitDist(), se simulan 2000 observaciones y se compara su distribución con la de los precios observados.

# Simulación de datos
simulacion <- rBCCG(
                n  = 2000,
                mu = 749.2743664,
                sigma = exp(-0.7520409),
                nu = 0.2530936
              )

ggplot() +
  geom_histogram(
    data = datos,
    aes(x = precio),
    alpha = 0.5, fill = "gray50") +
  geom_histogram(
    data = data.frame(simulacion),
    aes(x = simulacion),
    alpha = 0.3, fill = "firebrick") +
  labs(title = "Distribución real vs Distribución simulada BCCG") +
  theme_bw()

Distribuciones truncadas

Una distribución truncada es aquella en la que, el rango de posibles valores para las variables \(Y\), es un subconjunto del rango de otra distribución. No hay que confundir una distribución truncada con una censurada. En la primera, valores por encima o por debajo de un determinado valor no existen, mientras que en las segundas los valores sí existen pero no pueden ser observados.

La versión truncada de cualquier distribución presente en gamlss.family puede obtenerse gracias a la función gen.trun() del paquete gamlss.tr. La idea es sencilla: se parte una distribución ya existente en gamlss, se determinan los límites de truncado (izquierda, derecha o ambos) y automáticamente se generan nuevas funciones d, p, q y r con esas características.

En el siguiente ejemplo se parte de una distribución t-student (TF()) y se trunca en ambas colas con límites 0 y 100.

# Distribución T con parámetros:
# media: mu = 80,
# varianza: sigma = 20 
# grados de libertad: nu = 5
simulacion <- rTF(n = 2000, mu = 80, sigma = 20, nu = 5)

library(gamlss.tr)

# Misma distribución pero truncada entre 0 y 100.
gen.trun(
  par  = c(0,100),
  family = "TF",
  name = "0_100",
  type = "both"
)

## A truncated family of distributions from TF has been generated 
##  and saved under the names:  
##  dTF0_100 pTF0_100 qTF0_100 rTF0_100 TF0_100 
## The type of truncation is both 
##  and the truncation parameter is 0 100

simulacion_truncada <- rTF0_100(n = 2000, mu = 80, sigma = 20, nu = 5)

p1 <- ggplot(data = data.frame(simulacion)) +
      geom_histogram(
        aes(x = simulacion, y = after_stat(density)),
        color = "black",
        alpha = 0.3) +
      stat_function(
        fun   = function(.x){dTF(x = .x, mu = 80, sigma = 20, nu = 5)},
        color = "darkblue",
        size  = 1) +
      labs(title = "Distribución original") +
      theme_bw()

p2 <- ggplot(data = data.frame(simulacion_truncada)) +
      geom_histogram(
        aes(x = simulacion_truncada, y = after_stat(density)),
        color = "black",
        alpha = 0.3) +
      stat_function(
        fun   = function(.x){dTF0_100(x = .x, mu = 80, sigma = 20, nu = 5)},
        color = "darkblue",
        size  = 1) +
      labs(title = "Distribución truncada (0, 100)") +
      theme_bw()

ggarrange(p1, p2)

Distribuciones censuradas

Una distribución censurada es aquella en la que, el valor de las variables \(Y\) se desconoce cuando es superior o inferior a un determinado límite. La diferencia con las distribuciones truncadas es que, en estas últimas, los valores por debajo o arriba de un determinado valor no existen, mientras que en las distribuciones censuradas, los valores sí existen pero se desconocen.

La versión censurada de cualquier distribución presente en gamlss.family puede obtenerse gracias a la función gen.cens() del paquete gamlss.tr. La idea es sencilla: se parte una distribución ya existente en gamlss, se determinan los límites de censura (izquierda, derecha o intervalo) y automáticamente se generan nuevas funciones d, p, q y r con esas características.

Distribuciones mixtas

Las distribuciones mixtas son aquellas que se crean como resultado de combinar varias distribuciones. Este tipo de ajustes son útiles cuando la distribución de la variable respuesta es de tipo multimodal o para crear modelos con variables latentes.

Cualquier combinación de distribuciones gamlss.family puede ser ajustada empleando dos funciones del paquete gamlss.mx:

• gamlssMX(): cuando las distribuciones no tienen parámetros en común, pudiendo incluso ser distribuciones de distintas familias.

• gamlssNP(): cuando las distribuciones tienen parámetros en común.

En ambos casos, el modelo se ajusta empleando el algoritmo EM.

library(gamlss.mx)
data(enzyme)

ggplot(data = enzyme, aes(x = act)) +
  geom_histogram(color = "black", alpha = 0.3, bins = 50) +
  geom_rug() +
  labs(title = "Distribución actividad enzimática") +
  theme_bw()

modelo_mx <- gamlssMX(
                formula = act ~ 1,
                data    = enzyme,
                family  = RG,
                K       = 2,
                control = MX.control(plot = FALSE)
             )
modelo_mx

## 
## Mixing Family:  c("RG", "RG") 
## 
## Fitting method: EM algorithm 
## 
## Call:  gamlssMX(formula = act ~ 1, family = RG, K = 2, data = enzyme,  
##     control = MX.control(plot = FALSE)) 
## 
## Mu Coefficients for model: 1 
## (Intercept)  
##       1.127  
## Sigma Coefficients for model: 1 
## (Intercept)  
##      -1.091  
## Mu Coefficients for model: 2 
## (Intercept)  
##      0.1557  
## Sigma Coefficients for model: 2 
## (Intercept)  
##      -2.641  
## 
## Estimated probabilities: 0.376018 0.623982 
## 
## Degrees of Freedom for the fit: 5 Residual Deg. of Freedom   240 
## Global Deviance:     86.2916 
##             AIC:     96.2916 
##             SBC:     113.798

d_modelo_mx <- getpdfMX(modelo_mx)
d_modelo_mx(y = 2)

## [1] 0.0772

ggplot(data = enzyme) +
geom_histogram(
  aes(x = act, y = after_stat(density)),
  color = "black",
  alpha = 0.3) +
stat_function(
  fun   = function(.x){d_modelo_mx(y = .x)},
  color = "darkblue",
  size  = 1) +
labs(title = "Distribución estimada por modelo mixto \n de dos componentes") +
theme_bw()

modelo_mx <- gamlssMXfits(
                n       = 10,
                formula = act ~ 1,
                data    = enzyme,
                family  = RG,
                K       = 2,
                control = MX.control(plot = FALSE)
             )

data(glasses)
glasses <- glasses %>%
           mutate(
             sex = recode_factor(sex, "1" = "hombre", "2" = "mujer")
           )

ggplot(data = glasses, aes(x = ageread)) +
  geom_histogram(alpha = 0.3, color = "black") +
  labs(title = "Distribución edad uso de gafas") +
  theme_bw() +
  theme(legend.position = "none")

ggplot(data = glasses, aes(x = ageread, fill = sex, color = sex)) +
  geom_histogram(alpha = 0.3, bins = 20) +
  labs(title = "Distribución edad uso de gafas") +
  facet_wrap(facets = vars(sex), nrow = 2, scales = "free") +
  theme_bw() +
  theme(legend.position = "none")

set.seed(123)
modelo_mx_no <- gamlssMX(
                  formula = ageread ~ sex,
                  data    = glasses,
                  family  = NO,
                  K       = 2,
                  control = MX.control(plot = FALSE)
               )

modelo_mx_ga <- gamlssMX(
                  formula = ageread ~ sex,
                  data    = glasses,
                  family  = GA,
                  K       = 2,
                  control = MX.control(plot = FALSE)
               )

GAIC(
  modelo_mx_no,
  modelo_mx_ga
)

d_modelo_mx <- getpdfMX(modelo_mx_ga)

p1 <- ggplot(data = glasses %>% filter(sex == "hombre")) +
      geom_histogram(
        aes(x = ageread, y = after_stat(density)),
        color = "black",
        alpha = 0.3) +
      stat_function(
        fun   = function(.x){d_modelo_mx(y = .x)},
        color = "darkblue",
        size  = 1) +
      facet_wrap(vars(sex)) +
      theme_bw()

p2 <- ggplot(data = glasses %>% filter(sex == "mujer")) +
      geom_histogram(
        aes(x = ageread, y = after_stat(density)),
        color = "black",
        alpha = 0.3) +
      stat_function(
        fun   = function(.x){d_modelo_mx(y = .x)},
        color = "darkblue",
        size  = 1) +
      facet_wrap(vars(sex)) +
      theme_bw()

ggarrange(p1, p2)

# plotZAGA(mu = 3, sigma = 0.5, nu = 0.5, from = 0, to = 10, 
#          n = 101,  main=NULL)

x <- seq(0, 10, length.out = 500)
y <- dZAGA(x = x, mu = 3, sigma = 0.5, nu = 0.5)
ggplot(data = data.frame(x, y), aes(x,y)) +
  geom_point() +
  labs(title = "Zero-adjusted gamma distribution")+ 
  theme_bw()

Datos

El set de datos film90 del paquete gamlss.data contiene información sobre el número de proyecciones y la recaudación conseguida en su primera semana de estreno (ambas en escala logarítmica) para 4015 películas.

data("film90")
datos <- film90
datos <- film90 %>%
         dplyr::select(log_recaudacion = lborev1, log_proyecciones =  lboopen)

ggplot(data = datos, aes(x = log_proyecciones, y = log_recaudacion)) +
  geom_point(alpha = 0.3)  +
  labs(title = "Recaudación vs Número de proyecciones") +
  theme_bw()

A partir del gráfico se evidencia que:

• La relación entre la variable respuesta y el predictor no es lineal.

• La dispersión y forma de la distribución varía en los diferentes niveles de la variable respuesta.
  
  

P-Splines

La función pb() permite incorporar P-Splines (Penalized Smoothing Splines) a los predictores continuos del modelo. Esta función emplea el método de local maximum likelihood para seleccionar automáticamente los grados de libertad efectivos óptimos (flexibilidad).

modelo_ps <- gamlss(
              formula = log_recaudacion ~ pb(log_proyecciones),
              data    = datos,
              family  = NO,
              trace   = FALSE
            )

summary(modelo_ps)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = log_recaudacion ~ pb(log_proyecciones),  
##     family = NO, data = datos, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.347147   0.087053   26.96   <2e-16 ***
## pb(log_proyecciones) 0.928889   0.007149  129.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33120    0.01114   29.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  4031 
## Degrees of Freedom for the fit:  12.73672
##       Residual Deg. of Freedom:  4018.263 
##                       at cycle:  2 
##  
## Global Deviance:     14109.58 
##             AIC:     14135.05 
##             SBC:     14215.32 
## ******************************************************************

datos <- datos %>%
         mutate(prediccion_ps = fitted(modelo_ps))

ggplot(
  data = datos,
  aes(x = log_proyecciones, y = log_recaudacion)) +
geom_point(alpha = 0.2)  +
geom_line(
  aes(x = log_proyecciones, y = prediccion_ps, color = "P-splines"),
  size = 1) +
scale_color_manual(
  breaks = c("P-splines"),
  values = c("P-splines" = "red")) +
labs(title = "Recaudación vs Número de proyecciones") +
theme_bw() +
theme(legend.position = "bottom")

Cuando se incorpora en algún predictor de un modelo aditivo (GAM o GAMLSS) una función smooth, el coeficiente de regresión y su error estimado no se pueden interpretar de la forma convencional, ya que estos solo contemplan la aportación lineal, ignorando la no lineal. La mejor forma de interpretar el impacto de una función smooth no lineal es:

• Mediante gráficos de dependencia parcial, que muestran el impacto en las predicciones del modelo a medida que varía el valor de un predictor y se mantienen constantes los otros. getPEF()

• Reajustar el modelo excluyendo el predictor de interés y evaluar el impacto que tiene sobre el modelo. drop1()

getPEF(
  modelo_ps,
  term      = "log_proyecciones",
  parameter = "mu",
  plot      = TRUE
)

## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3)

drop1(object = modelo_ps, parallel = "multicore", ncpus = 4)

Monotonic P-Splines

La función pbm() incorpora una modificación de las P-splines que fuerza a la que curva resultante sea monótona (siempre ascendente o siempre descendente). Los argumentos son los mismo que los de pb() pero con el añadido de mono, que indica la dirección que debe tener la curva (“up” o “down”).

modelo_psm <- gamlss(
                    formula = log_recaudacion ~ pbm(log_proyecciones, mono = "up"),
                    data    = datos,
                    family  = NO,
                    trace   = FALSE
                  )

datos <- datos %>%
         mutate(prediccion_psm = fitted(modelo_psm))

ggplot(
  data = datos,
  aes(x = log_proyecciones, y = log_recaudacion)) +
geom_point(alpha = 0.2)  +
geom_line(
  aes(x = log_proyecciones, y = prediccion_ps, color = "P-splines"),
  size = 1) +
geom_line(
  aes(x = log_proyecciones, y = prediccion_psm, color = "Monotonic P-splines"),
  size = 1) +
scale_color_manual(
  breaks = c("P-splines", "Monotonic P-splines"),
  values = c("P-splines" = "red", "Monotonic P-splines" = "blue")) +
labs(title = "Recaudación vs Número de proyecciones") +
theme_bw() +
theme(legend.position = "bottom")

Cubic Splines

La función cs() permite incorporar Cubic Splines a los predictores continuos del modelo. A diferencia de pb(), en este caso sí es necesario especificar los grados de libertad.

# Modelo con cubic splines de 5 grados de libertad.
modelo_cs <- gamlss(
              formula = log_recaudacion ~ cs(log_proyecciones, df = 5),
              data    = datos,
              family  = NO,
              trace   = FALSE
            )

summary(modelo_cs)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = log_recaudacion ~ cs(log_proyecciones,  
##     df = 5), family = NO, data = datos, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  2.347158   0.087369   26.86   <2e-16 ***
## cs(log_proyecciones, df = 5) 0.928888   0.007175  129.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.33481    0.01114   30.06   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  4031 
## Degrees of Freedom for the fit:  8.000946
##       Residual Deg. of Freedom:  4022.999 
##                       at cycle:  2 
##  
## Global Deviance:     14138.75 
##             AIC:     14154.75 
##             SBC:     14205.17 
## ******************************************************************

datos <- datos %>%
         mutate(prediccion_cs = fitted(modelo_cs))

ggplot(
  data = datos,
  aes(x = log_proyecciones, y = log_recaudacion)) +
geom_point(alpha = 0.2)  +
geom_line(
  aes(x = log_proyecciones, y = prediccion_ps, color = "P-splines"),
  size = 1) +
geom_line(
  aes(x = log_proyecciones, y = prediccion_cs, color = "Cubic-splines"),
  size = 1) +
scale_color_manual(
  breaks = c("P-splines", "Cubic-splines"),
  values = c("P-splines" = "red", "Cubic-splines" = "blue")) +
labs(title = "Recaudación vs Número de proyecciones") +
theme_bw() +
theme(legend.position = "bottom")

Redes neuronales

Las redes neuronales (multiperceptrón) también pueden emplearse como funciones smooth (requiere del paquete gamlss.add).

library(gamlss.add)
modelo_nn <- gamlss(
              formula = log_recaudacion ~ nn(~log_proyecciones, size = 20, decay = 0.1),
              data    = datos,
              family  = NO,
              trace   = FALSE
            )

datos <- datos %>%
         mutate(prediccion_nn = fitted(modelo_nn))

ggplot(
  data = datos,
  aes(x = log_proyecciones, y = log_recaudacion)) +
geom_point(alpha = 0.2)  +
geom_line(
  aes(x = log_proyecciones, y = prediccion_ps, color = "P-splines"),
  size = 1) +
geom_line(
  aes(x = log_proyecciones, y = prediccion_cs, color = "Cubic-splines"),
  size = 1) +
geom_line(
  aes(x = log_proyecciones, y = prediccion_nn, color = "red-neuronal"),
  size = 0.7) +
scale_color_manual(
  breaks = c("P-splines", "Cubic-splines", "red-neuronal"),
  values = c("P-splines"="red", "Cubic-splines"="blue", "red-neuronal"="green")) +
labs(title = "Recaudación vs Número de proyecciones") +
theme_bw() +
theme(legend.position = "bottom")

Cando el modelo tiene un único predictor, resulta útil para su interpretación representar cómo varía la distribución estimada para la variable respuesta en función del predictor. La función plotSimpleGamlss() del paquete gamlss.util genera este tipo de gráficos.

modelo_gamlss <- gamlss(
                  formula = log_recaudacion ~ pb(log_proyecciones),
                  sigma.formula = log_recaudacion ~ pb(log_proyecciones),
                  data    = datos,
                  family  = NO,
                  trace   = FALSE
                )

library(gamlss.util)
library(colorspace)

plotSimpleGamlss(
  y = log_recaudacion,
  x = log_proyecciones,
  model = modelo_gamlss,
  data  = datos,
  x.val = seq(6, 16, 2),
  val   = 5,
  N     =  1000,
  ylim  = c(0,25),
  cols  = heat_hcl(100)
)

## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3) 
## new prediction 
## New way of prediction in pb()  (starting from GAMLSS version 5.0-3)

Error de los coeficientes

El paquete gamlss incorpora dos funciones para calcular los errores estándar de los parámetros estimados: vcov() y rvcov(). La primera es la utilizada por defecto en el summary() del modelo, y la segunda es una estimación robusta (robust or sandwich standard errors). El orden en el que aparecen los resultados de cada parámetro es \((\mu, \sigma, \nu, \tau)\), aunque no todos los modelos tienen que tener los 4.

modelo_gamlss <- gamlss(
                  formula = log_recaudacion ~ log_proyecciones,
                  sigma.formula = log_recaudacion ~ log_proyecciones,
                  data    = datos,
                  family  = NO,
                  trace   = FALSE
                )

vcov(modelo_gamlss, type = "se")

##      (Intercept) log_proyecciones      (Intercept) log_proyecciones 
##      0.106965570      0.007378298      0.054075872      0.004490221

rvcov(modelo_gamlss, type = "se")

##      (Intercept) log_proyecciones      (Intercept) log_proyecciones 
##      0.106793092      0.007090207      0.047121586      0.004005457

En cualquiera de los dos casos, para predictores que incluyan funciones smooth, su coeficiente solo contempla la contribución lineal, por lo tanto, no deben considerarse o interpretarse con mucha cautela.

Significancia

El p-value asociado a cada uno de los términos del modelo mostrado en el summary se calcula mediante el Wald test. Los propios autores del paquete recomiendan que, si este valor es de importancia para el analista, se emplee la función drop1() en su lugar. Esta función calcula, mediante un generalized likelihood ratio test (GLRT), la significancia del impacto que tiene excluir del modelo cada término. En ambos casos, las estimaciones solo aplican si el modelo no incluye términos con funciones smooth. En el caso de que sí los incluya, el resultado obtenido solo sirve como una aproximación.

modelo_gamlss <- gamlss(
                  formula       = log_recaudacion ~ log_proyecciones,
                  sigma.formula = log_recaudacion ~ log_proyecciones,
                  data    = datos,
                  family  = NO,
                  trace   = FALSE
                )

drop1(object = modelo_gamlss, parallel = "multicore", ncpus = 4)

Intervalos de confianza

EL paquete gamlss incorpora tres funciones distintas para calcular el intervalo de confianza de los parámetros estimados. Los resultados mostrados en el summary se corresponden con los calculados por la función confint(), que emplea intervalos de confianza basados en el Wald standard error. Las dos otras alternativas son las funciones prof.dev() y prof.term() que calculan intervalos basados en el profile likelihood.

modelo_gamlss <- gamlss(
                  formula       = log_recaudacion ~ log_proyecciones,
                  sigma.formula = log_recaudacion ~ log_proyecciones,
                  data    = datos,
                  family  = NO,
                  trace   = FALSE
                )

La función confint() requiere que se especifique sobre qué parámetro se quiere calcular el intervalo de confianza.

confint(
  modelo_gamlss,
  what   = c("mu"),
  level  = 0.95,
  robust = TRUE
)

##                      2.5 %    97.5 %
## (Intercept)      2.9031354 3.3217566
## log_proyecciones 0.8548881 0.8826812

confint(
  modelo_gamlss,
  what   = c("sigma"),
  level  = 0.95,
  robust = TRUE
)

##                       2.5 %     97.5 %
## (Intercept)       1.9732621  2.3918833
## log_proyecciones -0.1726126 -0.1448195

Predicción

Una vez que el modelo gamlss ha sido entrenado, se puede predecir el valor de cada uno de los parámetros para nuevos valores de los predictores. Dos funciones pueden emplearse con este fin:

• predict(): genera las predicciones para uno de los parámetros del modelo, especificado por el usuario.

• predictAll(): genera las predicciones para todos parámetros del modelo.

Ambas funciones permiten calcular tres tipos distintos de predicciones. Su comportamiento se controla con el argumento type:

• type="link": devuelve la predicción del parámetro en la escala de la función link. Por ejemplo, si la función link de la varianza del modelo es \(\log\), el valor devuelto es el logaritmo de la varianza.

• type="response": devuelve la predicción del parámetro en la misma escala de la variable respuesta.

• type="terms": muestra la contribución de cada término del modelo a la predicción.

modelo_gamlss <- gamlss(
                  formula = log_recaudacion ~ log_proyecciones,
                  sigma.formula = log_recaudacion ~ log_proyecciones,
                  data    = datos,
                  family  = NO,
                  trace   = FALSE
                )

Para la media \(\mu\) de la distribución normal, como su función link es la identidad, los resultados obtenidos empleando type="link" y type="response" son iguales.

predict(
  object  = modelo_gamlss,
  what    = "mu",
  type    = "link",
  newdata = data.frame(log_proyecciones = 11)
)

## [1] 12.66908

predict(
  object  = modelo_gamlss,
  what    = "mu",
  type    = "response",
  newdata = data.frame(log_proyecciones = 11)
)

## [1] 12.66908

predict(
  object  = modelo_gamlss,
  what    = "mu",
  type    = "terms",
  newdata = data.frame(log_proyecciones = 11)
)

##   log_proyecciones
## 1       -0.6808943
## attr(,"constant")
## [1] 13.34997

En el caso de la varianza, la función link empleada es \(\log\), por lo que el resultado si se indica type="response" equivale al exponente del resultado obtenido con type="link".

predict(
  object  = modelo_gamlss,
  what    = "sigma",
  type    = "link",
  newdata = data.frame(log_proyecciones = 11)
)

## [1] 0.4366962

predict(
  object  = modelo_gamlss,
  what    = "sigma",
  type    = "response",
  newdata = data.frame(log_proyecciones = 11)
)

## [1] 1.547586

La función predictAll() funciona del mismo modo, pero devuelve la predicción de todos los parámetros de la distribución a la vez.

predictAll(
  object  = modelo_gamlss,
  type    = "response",
  newdata = data.frame(log_proyecciones = 11)
)

## $mu
## [1] 12.66908
## 
## $sigma
## [1] 1.547586
## 
## attr(,"family")
## [1] "NO"     "Normal"

A la hora de generar predicciones, es muy importante tener en cuenta si el valor de los predictores está dentro del rango observado durante el entrenamiento del modelo. De no ser así, hay que ser muy cauto, ya que la relación aprendida entre la variable respuesta y los predictores podría no cumplirse fuera de las regiones observadas. Es buena práctica evitar las predicciones en regiones extrapoladas.

Métricas de ajuste

Tres métodos están disponibles en gamlss para cuantificar la capacidad predictiva de los modelos: Generalized Akaike information criterion (GAIC), validación simple y cross-validation.

data("rent")
datos <- rent
datos <- datos %>% select(R, Fl, A, H, loc)
# Se renombran las variables para que sean más explicativas
colnames(datos) <- c("precio", "metros", "anyo", "calefaccion", "situacion")

# Modelos candidatos
modelo_1 <- gamlss(formula = precio ~ metros, 
                   family = NO, data = datos, trace = FALSE)
modelo_2 <- gamlss(formula = precio ~ metros + anyo,
                   family = NO, data = datos, trace = FALSE)
modelo_3 <- gamlss(formula = precio ~ metros,
                   family = GA, data = datos, trace = FALSE)
modelo_4 <- gamlss(formula = precio ~ metros + anyo,
                   family = GA, data = datos, trace = FALSE)

GAIC(modelo_1, modelo_2, modelo_3, modelo_4, k = 3)

cv_modelo_1 <- gamlssCV(formula = precio ~ metros, family = NO, data = datos,
                        K.fold = 10, parallel = "multicore", ncpus = 4,
                        set.seed = 1)
cv_modelo_2 <- gamlssCV(formula = precio ~ metros + anyo, family = NO, data = datos,
                        K.fold = 10, parallel = "multicore", ncpus = 4,
                        set.seed = 1)
cv_modelo_3 <- gamlssCV(formula = precio ~ metros, family = GA, data = datos,
                        K.fold = 10, parallel = "multicore", ncpus = 4,
                        set.seed = 1)
cv_modelo_4 <- gamlssCV(formula = precio ~ metros + anyo, family = GA, data = datos,
                        K.fold = 10, parallel = "multicore", ncpus = 4,
                        set.seed = 1)

CV(cv_modelo_1, cv_modelo_2, cv_modelo_3, cv_modelo_4)

# Partición de los datos
set.seed(123)
id_train <- sample(x = 1:nrow(datos), size = nrow(datos)*0.8, replace = FALSE)
id_test  <- (1:nrow(datos))[-id_train]

# Modelos candidatos entrenados con la partición de entrenamiento
modelo_1 <- gamlss(formula = precio ~ metros, 
                   family = NO, data = datos[id_train, ], trace = FALSE)
modelo_2 <- gamlss(formula = precio ~ metros + anyo,
                   family = NO, data = datos[id_train, ], trace = FALSE)
modelo_3 <- gamlss(formula = precio ~ metros,
                   family = GA, data = datos[id_train, ], trace = FALSE)
modelo_4 <- gamlss(formula = precio ~ metros + anyo,
                   family = GA, data = datos[id_train, ], trace = FALSE)

# Calculo métrica de test
test_modelo_1 <- getTGD(modelo_1, newdata = datos[id_test, ])
test_modelo_2 <- getTGD(modelo_2, newdata = datos[id_test, ])
test_modelo_3 <- getTGD(modelo_3, newdata = datos[id_test, ])
test_modelo_4 <- getTGD(modelo_4, newdata = datos[id_test, ])
TGD(test_modelo_1, test_modelo_2, test_modelo_3, test_modelo_4)

Método de búsqueda

En el apartado anterior, se ha mostrado cómo comparar modelos asumiendo que se conoce cuáles son los candidatos a comparar. Cuando no se sabe qué predictores pueden ser los adecuados, conviene seguir una estrategia de búsqueda.

Una complejidad añadida de los modelos GAMLSS es que, cada parámetro del modelo, puede modelarse empleando distintos predictores. Por lo tanto, la búsqueda debe realizarse una vez por cada parámetro que tenga la distribución seleccionada. Las funciones stepGAICAll.A(), stepGAICAll.B() y stepTGDAll.A() ayudan en este proceso.

stepGAICAll.A()

modelo_nulo <- gamlss(precio~1, data = datos, family = GA, trace = FALSE)
mejor_modelo <- stepGAICAll.A(
                  object = modelo_nulo,
                  scope  = list(lower=~1, upper=~metros+anyo+calefaccion+situacion),
                  k      = log(nrow(datos)),
                  parallel = "multicore",
                  ncpus  = 4
                )

## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 28626.75 
##  precio ~ 1 
## 
##               Df   AIC
## + metros       1 28121
## + calefaccion  1 28422
## + situacion    2 28567
## + anyo         1 28593
## <none>           28627
## 
## Step:  AIC= 28121.18 
##  precio ~ metros 
## 
##               Df   AIC
## + calefaccion  1 27881
## + anyo         1 28018
## + situacion    2 28054
## <none>           28121
## 
## Step:  AIC= 27881.43 
##  precio ~ metros + calefaccion 
## 
##             Df   AIC
## + situacion  2 27833
## + anyo       1 27864
## <none>         27881
## 
## Step:  AIC= 27832.68 
##  precio ~ metros + calefaccion + situacion 
## 
##        Df   AIC
## + anyo  1 27818
## <none>    27833
## 
## Step:  AIC= 27817.69 
##  precio ~ metros + calefaccion + situacion + anyo 
## 
## --------------------------------------------------- 
## Distribution parameter:  sigma 
## Start:  AIC= 27817.69 
##  ~1 
## 
##               Df   AIC
## + anyo         1 27782
## + calefaccion  1 27808
## <none>           27818
## + metros       1 27821
## + situacion    2 27824
## 
## Step:  AIC= 27781.76 
##  ~anyo 
## 
##               Df   AIC
## <none>           27782
## + calefaccion  1 27786
## + metros       1 27788
## + situacion    2 27789
## --------------------------------------------------- 
## Distribution parameter:  mu 
## Start:  AIC= 27781.76 
##  precio ~ metros + calefaccion + situacion + anyo 
## 
##               Df   AIC
## <none>           27782
## - anyo         1 27807
## - situacion    2 27834
## - calefaccion  1 27902
## - metros       1 28383
## ---------------------------------------------------

mejor_modelo

## 
## Family:  c("GA", "Gamma") 
## Fitting method: RS() 
## 
## Call:  gamlss(formula = precio ~ metros + calefaccion + situacion +  
##     anyo, sigma.formula = ~anyo, family = GA, data = datos,      trace = FALSE) 
## 
## 
## Mu Coefficients:
##  (Intercept)        metros  calefaccion1    situacion2    situacion3  
##      1.92187       0.01082      -0.28960       0.20406       0.27559  
##         anyo  
##      0.00198  
## Sigma Coefficients:
## (Intercept)         anyo  
##    5.971699    -0.003575  
## 
##  Degrees of Freedom for the fit: 8 Residual Deg. of Freedom   1961 
## Global Deviance:     27721.1 
##             AIC:     27737.1 
##             SBC:     27781.8

El mejor modelo encontrado con stepGAICAll.A() emplea los predictores metros, calefaccion, situacion y anyo para modelar el parámetro \(\mu\), y únicamente el predictor anyo para modelar el parámetro \(\sigma\).

stepGAICAll.B()

La función stepGAICAll.B() sigue la misma estrategia que stepGAICAll.A() pero obligando a que todos los parámetros de la distribución sean modelados con son los mismos predictores.

modelo_nulo <- gamlss(precio~1, data = datos, family = GA, trace = FALSE)
mejor_modelo <- stepGAICAll.B(
                  object = modelo_nulo,
                  scope  = list(lower=~1, upper=~metros+anyo+calefaccion+situacion),
                  k      = log(nrow(datos)),
                  parallel = "multicore",
                  ncpus  = 4
                )

## Start:  AIC= 28626.75 
##  precio ~ 1 
## 
##               Df   AIC
## + metros       2 28121
## + calefaccion  2 28424
## + anyo         2 28576
## + situacion    4 28580
## <none>           28627
## 
## Step:  AIC= 28121.2 
##  precio ~ metros 
## 
##               Df   AIC
## + calefaccion  2 27871
## + anyo         2 27973
## + situacion    4 28060
## <none>           28121
## - metros       2 28627
## 
## Step:  AIC= 27870.78 
##  precio ~ metros + calefaccion 
## 
##               Df   AIC
## + situacion    4 27835
## + anyo         2 27841
## <none>           27871
## - calefaccion  2 28121
## - metros       2 28424
## 
## Step:  AIC= 27835.07 
##  precio ~ metros + calefaccion + situacion 
## 
##               Df   AIC
## + anyo         2 27800
## <none>           27835
## - situacion    4 27871
## - calefaccion  2 28060
## - metros       2 28397
## 
## Step:  AIC= 27799.49 
##  precio ~ metros + calefaccion + situacion + anyo 
## 
##               Df   AIC
## <none>           27800
## - anyo         2 27835
## - situacion    4 27841
## - calefaccion  2 27914
## - metros       2 28401

mejor_modelo

## 
## Family:  c("GA", "Gamma") 
## Fitting method: RS() 
## 
## Call:  gamlss(formula = precio ~ metros + calefaccion + situacion +  
##     anyo, sigma.formula = ~metros + calefaccion + situacion +  
##     anyo, family = GA, data = datos, trace = FALSE) 
## 
## Mu Coefficients:
##  (Intercept)        metros  calefaccion1    situacion2    situacion3  
##     1.911247      0.010903     -0.287947      0.201223      0.270280  
##         anyo  
##     0.001984  
## Sigma Coefficients:
##  (Intercept)        metros  calefaccion1    situacion2    situacion3  
##     4.934544      0.001269      0.075659     -0.106901     -0.157943  
##         anyo  
##    -0.003038  
## 
##  Degrees of Freedom for the fit: 12 Residual Deg. of Freedom   1957 
## Global Deviance:     27708.5 
##             AIC:     27732.5 
##             SBC:     27799.5

Ejemplo

El ajuste de modelos con gamboostLSS es muy similar al realizado con gamlss. Para mostrarlo, se emplea un set de datos india con el que se crea un modelo predictivo del índice de malnutrición infantil.

El índice stunted growth se define como la diferencia estandarizada entre el crecimiento de un niño y el crecimiento medio (mediana) de todos los niños de la población a la que pertenece. Está considerado como una de las principales consecuencias de la malnutrición de las madres durante el embarazo o de los niños durante los primeros meses de vida. Valores negativos indican un retraso en el crecimiento del niño.

Se crea un modelo para predecir el stunted growth en función de la edad e índice de masa corporal de la madre y de los niños durante los primeros meses de vida.

library(gamboostLSS)
data(india)
datos <- india

head(datos)

El paquete gamboostLSS tiene dos funciones para el ajuste de modelos: glmboostLSS() para modelos lineales y gamboostLSS() para todos los demás. En ambos casos, los argumentos principales son:

• formula: formula que define la variable respuesta y los predictores. Si solo se indica una fórmula, se emplea esa misma para estimar todos los parámetros de la distribución. También es posible pasar una lista en la que se especifica el modelo para cada uno de los parámetros.

• families: una lista con la familia de distribución para cada uno de los parámetros.

• boost-control(): parámetros para el algoritmo de boosting. De entre ellos destaca mstop que controla el ńumero de iteraciones de entrenamiento para cada parámetro, y nu que indica el ratio de aprendizaje (learning rate). Si solo se indica un valor de mstops se emplea ese número de iteraciones para ajustar todos los parámetros del modelo.

Aunque el paquete de gamboostLSS incorpora las principales distribuciones, se puede convertir cualquier distribución del paquete gamlss.dist con la función as.families(). Por ejemplo, para utilizar la función gamma se obtiene con as.families("GA").

En cuanto a los base learners disponibles para el ajuste por boosting, gamboostLSS emplea los disponibles en el paquete mboost:

• bbs(): para incorporar relaciones no lineales mediante penalized regression splines.
• bols(): para relaciones lineales.

Se procede a ajustar un modelo GAMLSS que prediga el índice stunting en función de la edad e índice de masa corporal del niño y de la madre. Se emplea una distribución normal para la variable respuesta, se modelan sus dos parámetros (\(\mu\) y \(\sigma\)) en función de los cuatro predictores, y se emplean penalized regression splines en todos los predictores para incorporar relaciones no lineales.

modelo <- gamboostLSS(
            formula = list(
                       mu = stunting ~ bbs(mage) + bbs(mbmi) + bbs(cage) + bbs(cbmi),
                       sigma = stunting ~ bbs(mage) + bbs(mbmi) + bbs(cage) + bbs(cbmi)
                      ),
            families = GaussianLSS(stabilization = "MAD"),
            data = india,
            control = boost_control(
                        trace = FALSE,
                        mstop = c(mu = 200, sigma = 100)
                      )
          )

Con el argumento mstop = c(mu = 200, sigma = 100) se está indicando que durante el ajuste, se vaya actualizando cada parámetro de forma alterna hasta alcanzar las 100 iteraciones A partir de ahí, el valor de sigma se mantiene constante y solo se sigue actualizando el de mu hasta alcanzar las 200 iteraciones.

summary(modelo)

## 
##   LSS Models fitted via Model-based Boosting
## 
## Call:
## gamboostLSS(formula = list(mu = stunting ~ bbs(mage) + bbs(mbmi) +     bbs(cage) + bbs(cbmi), sigma = stunting ~ bbs(mage) + bbs(mbmi) +     bbs(cage) + bbs(cbmi)), data = india, families = GaussianLSS(stabilization = "MAD"),     control = boost_control(trace = FALSE, mstop = c(mu = 200,         sigma = 100)))
## 
## Number of boosting iterations (mstop):  mu = 200, sigma = 100 
## Step size:  mu = 0.1, sigma = 0.1 
## 
## Families:
## 
##   Normal distribution: mu(id link) 
## 
## Loss function: -dnorm(x = y, mean = f, sd = sigma, log = TRUE) 
##  
## 
##   Normal distribution: sigma (log link) 
## 
## Loss function: -dnorm(x = y, mean = mu, sd = exp(f), log = TRUE) 
##  
## Selection frequencies:
## Parameter mu:
## bbs(cage) bbs(cbmi) bbs(mbmi) bbs(mage) 
##     0.460     0.235     0.160     0.145 
## Parameter sigma:
## bbs(mage) bbs(cage) bbs(cbmi) bbs(mbmi) 
##      0.40      0.21      0.21      0.18

Como en todo modelo ajustado empleando gradient boosting, el principal hiperparámetro es el número de iteraciones de entrenamiento mstop, ya que esto impacta directamente en el balance entre generalización y overfitting. En este tipo de modelos, no hay una forma clara de definir los grados de libertad, por lo que no puede emplearse como métrica de validación el GAIC. En su lugar, se recurre a validación cruzada (cross-validation) o bootstrapping como forma de cuantificar la bondad de ajuste de un modelo gamboostlss.

La optimización del hiperparámetro mstop en los modelos gamboostlss tiene una complejidad añadida comparado a los modelos de gradient boosting tradicionales. En este caso existe un valor de mstop distinto para cada parámetro de la distribución del modelo, por lo que es necesario una búsqueda combinada de todos ellos (grid multidimensional).

Todo este proceso puede hacerse con las funciones cvrisk() y make.grid(). En esta búsqueda es importante tener en cuenta la jerarquía de importancia en los parámetros del modelo. En términos generales, el parámetro \(\mu\) de cualquier distribución suele ser más crítico que \(\sigma\), por lo que vale la pena ser más exhaustivo en la búsqueda del valor de mstop asociado a este parámetro. Con este objetivo, la función make.grid() acepta el argumento dense_mu_grid para dar más granularidad en el grid de este parámetro.

Se realiza la búsqueda del mejor valor de mstop para cada parámetro del modelo anterior.

# Reducir length.out en caso de tiempo de ejecución excesivo
grid <- make.grid(
          max = c(mu = 3000, sigma = 500),
          min = 50,
          length.out = 5, 
          dense_mu_grid = FALSE
        )

head(grid)

ggplot(data = grid, aes(x = mu, y = sigma)) +
  geom_point() +
  labs(title = "Número de iteraciones de ajuste para mu y sigma") +
  theme_bw()

result_grid <- cvrisk(
                object = modelo,
                folds = cv(weights = model.weights(modelo), type = "kfold", B = 5),
                grid = grid, 
                mc.cores = 4
              )

## Starting cross-validation...
## [fold]   [current mstop]

result_grid

## 
##   Cross-validated risk
##   gamboostLSS(formula = list(mu = stunting ~ bbs(mage) + bbs(mbmi) +      bbs(cage) + bbs(cbmi), sigma = stunting ~ bbs(mage) + bbs(mbmi) +      bbs(cage) + bbs(cbmi)), data = india, families = GaussianLSS(stabilization = "MAD"),      control = boost_control(trace = FALSE, mstop = c(mu = 200,          sigma = 100))) 
## 
##    50,50   139,50   387,50  1078,50  3000,50    50,89   139,89   387,89 
## 1.817459 1.811000 1.810080 1.810340 1.810504 1.817789 1.811279 1.810322 
##  1078,89  3000,89   50,158  139,158  387,158 1078,158 3000,158   50,281 
## 1.810628 1.810804 1.817695 1.810981 1.809992 1.810315 1.810507 1.817456 
##  139,281  387,281 1078,281 3000,281   50,500  139,500  387,500 1078,500 
## 1.810700 1.809929 1.810268 1.810471 1.817464 1.810625 1.810187 1.810572 
## 3000,500 
## 1.810748 
## 
##   Optimal number of boosting iterations: 387 281

plot(result_grid)

Con la función mstop() se extrae el valor óptimo de mstop encontrado en la búsqueda.

mejores_hiperparametros <- mstop(result_grid)
mejores_hiperparametros

##    mu sigma 
##   387   281

Una vez identificados, se reestrena el modelo empleando los valores óptimos con la función mstop() aplicada a un modelo previamente definido.

mstop(modelo) <- mstop(result_grid)
modelo

## 
##   LSS Models fitted via Model-based Boosting
## 
## Call:
## gamboostLSS(formula = list(mu = stunting ~ bbs(mage) + bbs(mbmi) +     bbs(cage) + bbs(cbmi), sigma = stunting ~ bbs(mage) + bbs(mbmi) +     bbs(cage) + bbs(cbmi)), data = india, families = GaussianLSS(stabilization = "MAD"),     control = boost_control(trace = FALSE, mstop = c(mu = 200,         sigma = 100)))
## 
## Number of boosting iterations (mstop):  mu = 387, sigma = 281 
## Step size:  mu = 0.1, sigma = 0.1 
## 
## Families:
## 
##   Normal distribution: mu(id link) 
## 
## Loss function: -dnorm(x = y, mean = f, sd = sigma, log = TRUE) 
##  
## 
##   Normal distribution: sigma (log link) 
## 
## Loss function: -dnorm(x = y, mean = mu, sd = exp(f), log = TRUE) 
##