GAMLSS
2023-10-13
#1.- INTRODUCCIÓN
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…).
##1.1. - ¿QUÉ SON LOS MODELOS GAMLSS
Los modelos GAMLSS asumen que la variable respuesta tiene una función de densidad definida por hasta 4 parámetros (μ,σ,ν,τ) 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(yi|μi,σi,νi,τi) , la única condición es que el logaritmo de la función de densidad (logf(yi|μi,σi,νi,τi) ) y su primera derivada respecto a cada uno de los parámetros puedan calcularse.
1.2.- 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.
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 xi con la media de la variable respuesta g(μ) no es directa, sino que se hace a través de una función f(xi).
La función fi(xi) puede ser cualquier función que recibe el valor del predictor xi 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 μ puede ser modelada en función de otras variables (predictores) y cuya varianza σ se calcula mediante una constante de dispersión y una función ν(μ) . 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 (μ,σ,ν,τ).
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.
##1.3.- AJUSTES DE LOS MODELOS GAMLSS
##1.3.- aJUSTES 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(μ)=β0+f1(x1)+f2(x2)+…+fp(xp)
en términos generales, el algoritmo de ajuste por backfitting funciona de la siguiente forma: Se inicia con un valor para todos los términos fi del modelo, por ejemplo poniéndolos todos a cero.
Se estima el valor del primer término f1 empleando los datos de entrenamiento, con el objetivo de predecir lo mejor posible la variable respuesta y .
Se estima el valor del segundo término f2 empleando residuos y−f1(x1) .
Se repite el paso 3, ajustando cada término con los residuos del anterior.
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.
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).
##1.4.- 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. 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.
#2.- DATOS (EDA)
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.
# Librerías
# --------------------------
library(gamlss)## Warning: package 'gamlss' was built under R version 4.2.3## Loading required package: splines## Loading required package: gamlss.data##
## Attaching package: 'gamlss.data'## The following object is masked from 'package:datasets':
##
## sleep## Loading required package: gamlss.dist## Warning: package 'gamlss.dist' was built under R version 4.2.3## Loading required package: nlme## Loading required package: parallel## ********** GAMLSS Version 5.4-18 **********## For more on GAMLSS look at https://www.gamlss.com/## Type gamlssNews() to see new features/changes/bug fixes.library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'tibble' was built under R version 4.2.3## Warning: package 'tidyr' was built under R version 4.2.3## Warning: package 'readr' was built under R version 4.2.3## Warning: package 'purrr' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## Warning: package 'forcats' was built under R version 4.2.3## Warning: package 'lubridate' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::collapse() masks nlme::collapse()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(ggpubr)## Warning: package 'ggpubr' was built under R version 4.2.3library(skimr)## Warning: package 'skimr' was built under R version 4.2.3data("rent")datos <- rentdatos <- datos %>% select(R, Fl, A, H, loc)colnames(datos) <- c("precio", "metros", "anyo", "calefaccion", "situacion")str(datos)## 'data.frame': 1969 obs. of 5 variables:
## $ precio : num 693 422 737 732 1295 ...
## $ metros : num 50 54 70 50 55 59 46 94 93 65 ...
## $ anyo : num 1972 1972 1972 1972 1893 ...
## $ calefaccion: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 1 1 1 ...
## $ situacion : Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 2 2 2 2 ...head(datos)## precio metros anyo calefaccion situacion
## 1 693.3 50 1972 0 2
## 2 422.0 54 1972 0 2
## 3 736.6 70 1972 0 2
## 4 732.2 50 1972 0 2
## 5 1295.1 55 1893 0 2
## 6 1195.9 59 1893 0 2summary(datos)## precio metros anyo calefaccion situacion
## Min. : 101.7 Min. : 30.00 Min. :1890 0:1580 1: 172
## 1st Qu.: 544.2 1st Qu.: 52.00 1st Qu.:1934 1: 389 2:1247
## Median : 737.8 Median : 67.00 Median :1957 3: 550
## Mean : 811.9 Mean : 67.73 Mean :1948
## 3rd Qu.:1022.0 3rd Qu.: 82.00 3rd Qu.:1972
## Max. :3000.0 Max. :120.00 Max. :1988skim(datos)| Name | datos |
| Number of rows | 1969 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| factor | 2 |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| calefaccion | 0 | 1 | FALSE | 2 | 0: 1580, 1: 389 |
| situacion | 0 | 1 | FALSE | 3 | 2: 1247, 3: 550, 1: 172 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| precio | 0 | 1 | 811.88 | 379.00 | 101.7 | 544.2 | 737.8 | 1022 | 3000 | ▇▇▂▁▁ |
| metros | 0 | 1 | 67.73 | 20.86 | 30.0 | 52.0 | 67.0 | 82 | 120 | ▆▇▇▅▂ |
| anyo | 0 | 1 | 1948.48 | 29.02 | 1890.0 | 1934.0 | 1957.0 | 1972 | 1988 | ▃▁▃▇▇ |
Representación gráfica
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 calefaccion") +
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 situacion") +
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)
)
# Representar la variable respuesta (dependiente)
# ---------------------------------------------------
ggplot(data = datos, aes(x = precio)) +
geom_density(alpha = 0.5, fill = "blue") +
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.
#3.- MODELOS
##3.1.- MODELO LINEAL
# 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
## ******************************************************************lineal <- lm(precio ~ metros + anyo + calefaccion + situacion, data = datos)
summary(lineal)##
## Call:
## lm(formula = precio ~ metros + anyo + calefaccion + situacion,
## data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -973.60 -201.21 -36.87 177.51 1777.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2775.0388 527.6476 -5.259 1.60e-07 ***
## metros 8.8394 0.3391 26.068 < 2e-16 ***
## anyo 1.4808 0.2677 5.532 3.60e-08 ***
## calefaccion1 -204.7596 19.4080 -10.550 < 2e-16 ***
## situacion2 134.0523 25.1727 5.325 1.12e-07 ***
## situacion3 209.5815 27.1632 7.716 1.90e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 308.9 on 1963 degrees of freedom
## Multiple R-squared: 0.3372, Adjusted R-squared: 0.3355
## F-statistic: 199.7 on 5 and 1963 DF, p-value: < 2.2e-16La 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"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)## Warning in wp(modelo_lm, ylim.all = 1): Some points are missed out
## increase the y limits using ylim.all
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.
3.2.- MODELO GLM
La distribución gamma tiene dos parámetros que tienen que ser estimados: media (μ ) y escala (σ). Por defecto, el paquete gamlss emplea para ambos la función log como link, ya que los dos parámetros tienen un dominio (0, ∞).
η1=g1(μ)=log(μ)
η2=g1(σ)=log(σ)
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”
# 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 diponibles 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.5681289 5.043 5.01e-07 ***
## metros 0.0106232 0.0004128 25.735 < 2e-16 ***
## anyo 0.0015100 0.0002886 5.232 1.85e-07 ***
## calefaccion1 -0.3000745 0.0231153 -12.982 < 2e-16 ***
## situacion2 0.1907641 0.0305203 6.250 5.01e-10 ***
## situacion3 0.2640828 0.0329197 8.022 1.77e-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
## ******************************************************************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
## ******************************************************************Vemos entonces que el nuevo modelo funciona mejor que el modelo anterior.
# Worm plot de los residuos
wp(modelo_glm, ylim.all = 0.5)## Warning in wp(modelo_glm, ylim.all = 0.5): Some points are missed out
## increase the y limits using ylim.all
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)## df AIC
## modelo_glm 7 27778.59
## modelo_lm 7 28173.00El 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.
##3.3.- Modelo GAM
# 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.5692315 5.420 6.70e-08 ***
## pb(metros) 0.0103084 0.0004031 25.573 < 2e-16 ***
## pb(anyo) 0.0014062 0.0002893 4.861 1.26e-06 ***
## calefaccion1 -0.3008111 0.0225869 -13.318 < 2e-16 ***
## situacion2 0.1886692 0.0299295 6.304 3.58e-10 ***
## situacion3 0.2719856 0.0322862 8.424 < 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
## ******************************************************************Vemos que son estadísticos significativos, otra vez.
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.
# Esta función puede tardar si hay muchos predictores no lineales o muchos datos.
drop1(modelo_gam, parallel = "multicore", ncpus = 4)## Single term deletions for
## mu
##
## Model:
## precio ~ pb(metros) + pb(anyo) + calefaccion + situacion
## Df AIC LRT Pr(Chi)
## <none> 27706
## pb(metros) 1.4680 28261 558.59 < 2.2e-16 ***
## pb(anyo) 4.3149 27798 101.14 < 2.2e-16 ***
## calefaccion 1.8445 27862 160.39 < 2.2e-16 ***
## situacion 2.0346 27770 68.02 1.825e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Otra 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 (en este caso, la variable respuesta es el PRECIO de la vivienda). 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 η=log(μ) (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)## Warning in wp(modelo_gam, ylim.all = 0.5): Some points are missed out
## increase the y limits using ylim.all
GAIC(modelo_lm, modelo_glm,modelo_gam)## df AIC
## modelo_gam 11.21547 27705.65
## modelo_glm 7.00000 27778.59
## modelo_lm 7.00000 28173.00El modelo GAM ha conseguido reducir todavía más el valor AIC, y la distribución de sus residuos ha mejorado ligeramente.
##3.4.- GENERALIZAR EL GAM
Existen 3 formas distintas de parametrizar la distribución gamma, una de ellas es mediante la media μ y un parámetro de escala θ (llamado σ en el paquete gamlss). Hasta el momento, se ha modelado únicamente μ en función de los predictores, asumiendo que la escala es constante y puede ser estimada a partir de la media.
#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.5823154 4.953 7.92e-07 ***
## pb(metros) 0.0105402 0.0003620 29.117 < 2e-16 ***
## pb(anyo) 0.0014977 0.0002955 5.068 4.40e-07 ***
## calefaccion1 -0.2918923 0.0242789 -12.022 < 2e-16 ***
## situacion2 0.1938688 0.0323320 5.996 2.40e-09 ***
## situacion3 0.2734326 0.0339177 8.062 1.30e-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.8645543 6.850 9.84e-12 ***
## pb(metros) 0.0019227 0.0007727 2.488 0.0129 *
## pb(anyo) -0.0035795 0.0004372 -8.187 4.78e-16 ***
## calefaccion1 0.0659289 0.0415081 1.588 0.1124
## situacion2 -0.1166325 0.0565230 -2.063 0.0392 *
## situacion3 -0.1702147 0.0608449 -2.798 0.0052 **
## ---
## 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
## ******************************************************************term.plot(modelo_gamlss, parameter = "sigma", pages = 1, ask = FALSE, rug = TRUE)
drop1(modelo_gamlss, parameter = "sigma", parallel = "multicore", ncpus = 4)## Single term deletions for
## sigma
##
## Model:
## ~pb(metros) + pb(anyo) + calefaccion + situacion
## Df AIC LRT Pr(Chi)
## <none> 27615
## pb(metros) 4.02694 27631 24.683 5.997e-05 ***
## pb(anyo) 3.87807 27659 52.167 1.067e-10 ***
## calefaccion 0.88335 27615 1.866 0.14788
## situacion 2.03694 27619 8.036 0.01872 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Worm plot de los residuos
wp(modelo_gamlss, ylim.all=0.5)
Vemos que hay mucha mayor concentración entre las dos lóineas
discontinuas.
plot(modelo_gamlss)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.001490387
## variance = 1.001031
## coef. of skewness = -0.138024
## coef. of kurtosis = 3.075731
## Filliben correlation coefficient = 0.9991937
## ******************************************************************GAIC(
modelo_lm,
modelo_glm,
modelo_gam,
modelo_gamlss
)## df AIC
## modelo_gamlss 22.25035 27614.78
## modelo_gam 11.21547 27705.65
## modelo_glm 7.00000 27778.59
## modelo_lm 7.00000 28173.00Acorde 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 μ y otro de escala σ. 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 μ y σ . 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 η y ν , pueden ser modelados también en función de los predictores. Es aquí donde reside el potencial de estos modelos.
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-Cole 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.5875497 3.453 0.000567 ***
## pb(metros) 0.0104079 0.0003753 27.732 < 2e-16 ***
## pb(anyo) 0.0019277 0.0002974 6.481 1.15e-10 ***
## calefaccion1 -0.3213906 0.0269135 -11.942 < 2e-16 ***
## situacion2 0.1853274 0.0345355 5.366 9.00e-08 ***
## situacion3 0.2742129 0.0361897 7.577 5.44e-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.8654297 7.688 2.36e-14 ***
## pb(metros) 0.0018721 0.0008063 2.322 0.0203 *
## pb(anyo) -0.0039673 0.0004381 -9.055 < 2e-16 ***
## calefaccion1 0.0819273 0.0430278 1.904 0.0571 .
## situacion2 -0.0851622 0.0624983 -1.363 0.1732
## situacion3 -0.1410790 0.0672213 -2.099 0.0360 *
## ---
## 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.223936 -0.980 0.3271
## pb(metros) -0.000983 0.002315 -0.425 0.6712
## pb(anyo) 0.001989 0.001626 1.224 0.2213
## calefaccion1 -0.286663 0.123453 -2.322 0.0203 *
## situacion2 -0.171143 0.182756 -0.936 0.3492
## situacion3 -0.084593 0.201950 -0.419 0.6754
## ---
## 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
## ******************************************************************Vemos aquí que sólo es significativo calefaccion1.
term.plot(modelo_gamlss_2, parameter = "nu", pages =1, ask = FALSE, rug = TRUE)
drop1(modelo_gamlss_2, parameter = "nu", parallel = "multicore", ncpus = 4)## Single term deletions for
## nu
##
## Model:
## ~pb(metros) + pb(anyo) + calefaccion + situacion
## Df AIC LRT Pr(Chi)
## <none> 27608
## pb(metros) 1.0234 27606 0.2085 0.65746
## pb(anyo) 1.7773 27608 3.6102 0.13685
## calefaccion 1.1919 27612 5.7613 0.02197 *
## situacion 1.8542 27605 0.6695 0.68075
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El estudio de los coeficientes apunta a que únicamente el predictor calefacción contribuye a modelar el parámetro ν.
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
)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 = ~calefaccion,
## family = BCCGo, data = datos, trace = FALSE)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.4024759 0.5601346 4.289 1.88e-05 ***
## pb(metros) 0.0105022 0.0003577 29.364 < 2e-16 ***
## pb(anyo) 0.0017267 0.0002837 6.086 1.39e-09 ***
## calefaccion1 -0.3279448 0.0265806 -12.338 < 2e-16 ***
## situacion2 0.2009112 0.0319161 6.295 3.79e-10 ***
## situacion3 0.2834923 0.0335282 8.455 < 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) 6.3424228 0.8739134 7.257 5.67e-13 ***
## pb(metros) 0.0017822 0.0007934 2.246 0.0248 *
## pb(anyo) -0.0037936 0.0004413 -8.597 < 2e-16 ***
## calefaccion1 0.0834511 0.0429909 1.941 0.0524 .
## situacion2 -0.1085064 0.0592138 -1.832 0.0670 .
## situacion3 -0.1589369 0.0639265 -2.486 0.0130 *
## ---
## 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) 0.51093 0.05645 9.051 < 2e-16 ***
## calefaccion1 -0.34599 0.11303 -3.061 0.00224 **
## ---
## 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: 23.77443
## Residual Deg. of Freedom: 1945.226
## at cycle: 6
##
## Global Deviance: 27556.05
## AIC: 27603.6
## SBC: 27736.38
## ******************************************************************plot(modelo_gamlss_2)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.001976413
## variance = 1.000507
## coef. of skewness = -0.01539539
## coef. of kurtosis = 2.910624
## Filliben correlation coefficient = 0.9997421
## ******************************************************************wp(modelo_gamlss_2)
Vemos prácticamente todo centrado o colocado en el centro. Todo
concentrado en el cero.
GAIC(
modelo_lm,
modelo_glm,
modelo_gam,
modelo_gamlss,
modelo_gamlss_2
)## df AIC
## modelo_gamlss_2 23.77443 27603.60
## modelo_gamlss 22.25035 27614.78
## modelo_gam 11.21547 27705.65
## modelo_glm 7.00000 27778.59
## modelo_lm 7.00000 28173.00A lo largo de este ejemplo, se han mostrado las ventajas derivadas de la flexbilidad 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.
#4.- DISTRIBUCIONES
Los modelos GAMLSS permiten ajustar prácticamente cualquier distribución paramétrica. Esto supone una gran libertad para crear modelos de regresión, pero también más responsabilidad a la hora de decidir qué distribución emplear. Afortunadamente, el paquete gamlss y sus extensiones incorporan algunas estrategias para, junto con el conocimiento del analista, ayudar a identificar la mejor distribució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 (−∞,∞)
- 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,∞]
- 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)
library(gamlss)
distribuciones <- fitDist(
datos$precio,
K = 2, # penalización similar al AIC,
type = "realplus",
trace = F,
try.gamlss = T
)##
|
| | 0%
|
|=== | 4%
|
|====== | 9%
|
|========= | 13%
|
|============ | 17%
|
|=============== | 22%
|
|================== | 26%
|
|===================== | 30%
|
|======================== | 35%
|
|=========================== | 39%
|
|============================== | 43%## Warning in MLE(ll2, start = list(eta.mu = eta.mu, eta.sigma = eta.sigma), :
## possible convergence problem: optim gave code=1 false convergence (8)##
|
|================================= | 48%## Warning in MLE(ll2, start = list(eta.mu = eta.mu, eta.sigma = eta.sigma), :
## possible convergence problem: optim gave code=1 false convergence (8)##
|
|===================================== | 52%## Warning in MLE(ll2, start = list(eta.mu = eta.mu, eta.sigma = eta.sigma), :
## possible convergence problem: optim gave code=1 false convergence (8)##
|
|======================================== | 57%
|
|=========================================== | 61%
|
|============================================== | 65%Error in solve.default(oout$hessian) :
## Lapack routine dgesv: system is exactly singular: U[3,3] = 0
##
|
|================================================= | 70%
|
|==================================================== | 74%
|
|======================================================= | 78%## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function## Warning in body(fun): argument is not a function## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function## Warning in body(fun): argument is not a function## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function## Warning in body(fun): argument is not a function## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function## Warning in body(fun): argument is not a function##
|
|========================================================== | 83%
|
|============================================================= | 87%Error in solve.default(oout$hessian) :
## Lapack routine dgesv: system is exactly singular: U[4,4] = 0
##
|
|================================================================ | 91%
|
|=================================================================== | 96%
|
|======================================================================| 100%distribuciones$fits %>%
enframe(name = "Dist", value = "GAIC") %>%
arrange(GAIC)## # A tibble: 23 × 2
## Dist GAIC
## <chr> <dbl>
## 1 BCCG 28614.
## 2 BCCGo 28614.
## 3 GG 28614.
## 4 BCPEo 28615.
## 5 BCPE 28615.
## 6 GIG 28616.
## 7 BCTo 28616.
## 8 GA 28616.
## 9 BCT 28616.
## 10 GB2 28616.
## # ℹ 13 more rowssummary(distribuciones)## *******************************************************************
## Family: c("BCCG", "Box-Cox-Cole-Green")
##
## Call: gamlssML(formula = y, family = DIST[i], K = 2)
##
## 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.4plot(distribuciones)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 6.206975e-08
## variance = 1.000508
## coef. of skewness = 0.001736345
## coef. of kurtosis = 3.034594
## Filliben correlation coefficient = 0.9992592
## ******************************************************************wp(distribuciones)
#** 5.- SIMULACIÓN**
BCCG()##
## GAMLSS Family: BCCG Box-Cox-Cole-Green
## Link function for mu : identity
## Link function for sigma: log
## Link function for nu : identityset.seed(1234)
# 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] 405.3305# Simulación de Montecarlo
#------------------------------------
simulacion <- rBCCG(
n = 1000,
mu = 749.2743664,
sigma = exp(-0.7520409),
nu = 0.2530936
)ggplot()+
geom_histogram(
data = datos,
aes(x=precio),
alpha = 0.5, fill = "blue")+
geom_histogram(
data = data.frame(simulacion),
aes(x=simulacion),
alpha = 0.3, fill = "red")+
labs(title = "Real vs simulación")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
theme_bw()+
theme(legend.positon = "bottom")## List of 98
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 11
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ title : NULL
## $ aspect.ratio : NULL
## $ axis.title : NULL
## $ axis.title.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.75points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.75points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.bottom : NULL
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.75points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y.left : NULL
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.75points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.2points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.2points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.bottom : NULL
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 1
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.2points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y.left : NULL
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.2points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.x : NULL
## $ axis.ticks.x.top : NULL
## $ axis.ticks.x.bottom : NULL
## $ axis.ticks.y : NULL
## $ axis.ticks.y.left : NULL
## $ axis.ticks.y.right : NULL
## $ axis.ticks.length : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom: NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.x.top : NULL
## $ axis.line.x.bottom : NULL
## $ axis.line.y : NULL
## $ axis.line.y.left : NULL
## $ axis.line.y.right : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ legend.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.key.size : 'simpleUnit' num 1.2lines
## ..- attr(*, "unit")= int 3
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.align : NULL
## $ legend.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.title.align : NULL
## $ legend.position : chr "right"
## $ legend.direction : NULL
## $ legend.justification : chr "center"
## $ legend.box : NULL
## $ legend.box.just : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "unit")= int 1
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.box.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ panel.background :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.border :List of 5
## ..$ fill : logi NA
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.spacing : 'simpleUnit' num 5.5points
## ..- attr(*, "unit")= int 8
## $ panel.spacing.x : NULL
## $ panel.spacing.y : NULL
## $ panel.grid :List of 6
## ..$ colour : chr "grey92"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major : NULL
## $ panel.grid.minor :List of 6
## ..$ colour : NULL
## ..$ linewidth : 'rel' num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major.x : NULL
## $ panel.grid.major.y : NULL
## $ panel.grid.minor.x : NULL
## $ panel.grid.minor.y : NULL
## $ panel.ontop : logi FALSE
## $ plot.background :List of 5
## ..$ fill : NULL
## ..$ colour : chr "white"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ plot.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.title.position : chr "panel"
## $ plot.subtitle :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : num 1
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 5.5points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption.position : chr "panel"
## $ plot.tag :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag.position : chr "topleft"
## $ plot.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ strip.background :List of 5
## ..$ fill : chr "grey85"
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ strip.background.x : NULL
## $ strip.background.y : NULL
## $ strip.clip : chr "inherit"
## $ strip.placement : chr "inside"
## $ strip.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey10"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 4.4points 4.4points 4.4points 4.4points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.x : NULL
## $ strip.text.x.bottom : NULL
## $ strip.text.x.top : NULL
## $ strip.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.y.left :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.y.right : NULL
## $ strip.switch.pad.grid : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ strip.switch.pad.wrap : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ legend.positon : chr "bottom"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUE##6.- 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=1000, mu = 80, sigma = 20, nu = 5)library(gamlss.tr)## Warning: package 'gamlss.tr' was built under R version 4.2.3gen.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 100simulacion_truncada <- rTF0_100(n=1000, mu = 80, sigma = 20, nu = 5)p1 <- ggplot(data = data.frame(simulacion)) +
geom_histogram(
aes(x = simulacion, y = after_stat(density)),
color = "blue",
alpha = 0.3) +
stat_function(
fun = function(.x){dTF(x=.x, mu = 80, sigma = 20, nu=5)},
color = "red",
size = 1) +
labs(title = "Real") + theme_bw()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p2 <- ggplot(data = data.frame(simulacion_truncada)) +
geom_histogram(
aes(x = simulacion_truncada, y = after_stat(density)),
color = "blue",
alpha = 0.3) +
stat_function(
fun = function(.x){dTF0_100(x=.x, mu = 80, sigma = 20, nu=5)},
color = "red",
size = 1) +
labs(title = "Truncada") + theme_bw()library(ggpubr)
ggarrange(p1, p2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
##7.- DISTRIBUCIONES CENSURADAS
Una distribución censurada es aquélla 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.
##8.- DISTRIBUCIONES MIXTAS
Las distribuciones mixtas son aquéllas 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.
#9.- DISTRIBUCIONES BIMODALES
El set de datos enzyme del paquete gamlss.mx contiene información sobre la actividad en sangre de una enzima en 245 pacientes.
library(gamlss.mx)## Warning: package 'gamlss.mx' was built under R version 4.2.3## Loading required package: nnetlibrary(nnet)data("enzyme")library(ggplot2)
ggplot(data = enzyme, aes(x=act))+
geom_histogram(color = "blue", alpha = 0.3, bins = 50)+
geom_rug()+
labs(title = "Dist")+
theme_bw()+
theme(legend.position = "bottom")
La distribución es bimodal, lo que apunta a que es resultado de combinar
dos distribuciones distintas. Se ajusta un modelo mixto con dos
componentes (k=2), en el que cada componente es una
distribución Reverse Gumbel (family = RG).
modelo_mx <- gamlssMX(
formula = act ~ 1,
data = enzyme,
family = RG,
k = 2,
control = MX.control(plot = T)
)
modelo_mx##
## Mixing Family: c("RG", "RG")
##
## Fitting method: EM algorithm
##
## Call: gamlssMX(formula = act ~ 1, family = RG, data = enzyme,
## control = MX.control(plot = T), k = 2)
##
## Mu Coefficients for model: 1
## (Intercept)
## 0.1557
## Sigma Coefficients for model: 1
## (Intercept)
## -2.641
## Mu Coefficients for model: 2
## (Intercept)
## 1.127
## Sigma Coefficients for model: 2
## (Intercept)
## -1.091
##
## Estimated probabilities: 0.6239804 0.3760196
##
## Degrees of Freedom for the fit: 5 Residual Deg. of Freedom 240
## Global Deviance: 86.2916
## AIC: 96.2916
## SBC: 113.798Como el modelo resultante es una combinación de distribuciones, para obtener la densidad estimada de una nueva observación, hay que calcular primero la densidad de cada distribución por separado. Con la función getpdfMX() se pueden crear funciones auxiliares que automatizan el proceso.
d_modelo_mx <- getpdfMX(modelo_mx)
d_modelo_mx(y = 2)## [1] 0.07720399ggplot(data=enzyme)+
geom_histogram(
aes(act, y = after_stat(density)),
color = "blue",
alpha = 0.3) +
stat_function(
fun = function(.x){d_modelo_mx(y=.x)},
color = "red",
size = 1) +
labs(title = "Distribución") +
theme_bw()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
El ajuste por EM se caracteriza por correr el riesgo de caer en mínimos
locales y no llegar a la solución óptima. Para evitar este problema, se
aconseja ajustar el modelo varias veces, en cada una partiendo de unos
valores iniciales distintos. La función gamlssMXfits() automatiza este
proceso ajustando el modelo n veces y devolviendo al final el que mejor
resultado ha conseguido (menor global deviance).
modelo_mx <- gamlssMXfits(
n = 10,
formula = act ~ 1,
data = enzyme,
family = RG,
K = 2,
control = MX.control(plot = T)
)
## model= 1
## model= 2
## model= 3
## model= 4
## model= 5
## model= 6
## model= 7
## model= 8
## model= 9
## model= 10#10.- MODELOS CON VARIABLES LATENTES
Los modelos de variables latentes se emplean en escenarios en los que se asume que la población objetivo está dividida en varias clases, pero se desconocen cuáles son. El objetivo de estos modelos es estimar la probabilidad que tiene cada observación de pertenecer a cada una de las clases. En la mayoría de casos reales, el número total de las clases se desconoce, por lo que tiene que ser estimado a partir de los datos.
El set de datos glasses del paquete gamlss.data contiene información sobre la edad a la que empezaron a utilizar gafas 1016 personas. Se dispone también del sexo de cada persona.
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 = "bottom")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
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")
Se ajustan dos modelos mixtos, uno con componentes noramles y otro con
componentes gamma, en los que la variable respuesta es ageread y se
incorpora a la variable sex como predictor.
set.seed(12345)
modelo_mx_no <- gamlssMX(
formula = ageread ~ sex,
data = glasses,
family = NO,
K = 2,
control = MX.control(plot =T)
)
modelo_mx_ga <- gamlssMX(
formula = ageread ~ sex,
data = glasses,
family = GA,
K = 2,
control = MX.control(plot = T)
)
GAIC(modelo_mx_no, modelo_mx_ga)## df AIC
## modelo_mx_ga 7 7922.612
## modelo_mx_no 7 7945.059d_modelo_mx <- getpdfMX(modelo_mx_ga)
p1 <- ggplot(data = glasses %>% filter(sex == "hombre")) +
geom_histogram(
aes(x = ageread, y = after_stat(density)),
color = "blue",
alpha = 0.3) +
stat_function(
fun = function(.x){d_modelo_mx(y = .x)},
color = "red",
size = 0.5) +
facet_wrap(vars(sex)) +
theme_bw()
p2 <- ggplot(data = glasses %>% filter(sex == "mujer")) +
geom_histogram(
aes(x = ageread, y = after_stat(density)),
color = "blue",
alpha = 0.3) +
stat_function(
fun = function(.x){d_modelo_mx(y = .x)},
color = "red",
size = 0.5) +
facet_wrap(vars(sex)) +
theme_bw()
ggarrange(p1, p2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Dentro del objeto gamlssMX se almacenan información del modelo final,
destacan: los modelos individuales que forman parte del modelo mixto
(\(models) y la probabilidad de cada
observación de pertenecer a cada uno de las componentes
(\)post.prob).
#11.- Zero Inflated y Zero - Adjusted
Las distribuciones zero-inflated y zero-adjusted son casos especiales de distribuciones mixtas en las que uno de los componentes es exactamente 0 con probabilidad 1 y el otro componente es una distribución continua. Un ejemplo de variable que tiene este comportamiento es la precipitación. Para los días que no llueve, la cantidad de precipitación es exactamente 0; sin embargo, cuando llueve, la precipitación acumulada sigue una distribución continua. El paquete gamlss ya tiene predefinidas algunas funciones de este tipo como ZAGA() y ZAIG()
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()
#12.- Funciones no lineales
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()
#13.- 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 éstos sólo 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
drop1(object = modelo_ps, parallel = "multicore", ncpus = 4)## Single term deletions for
## mu
##
## Model:
## log_recaudacion ~ pb(log_proyecciones)
## Df AIC LRT Pr(Chi)
## <none> 14135
## pb(log_proyecciones) 10.737 20956 6842.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nos indica que es significativo.
#** 14.- MOTONIC - P-Splines**
La función pbm() incorpora una modificación de las P-splines que fuerza a que la curva resultante sea monótona (siempre ascendente o siempre descendente). Los argumentos son los mismos 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
)
summary(modelo_psm)## ******************************************************************
## Family: c("NO", "Normal")
##
## Call: gamlss(formula = log_recaudacion ~ pbm(log_proyecciones,
## mono = "up"), family = NO, data = datos, trace = FALSE)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.29293 0.02194 605.8 <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.33163 0.01114 29.78 <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: 11.90142
## Residual Deg. of Freedom: 4019.099
## at cycle: 2
##
## Global Deviance: 14113.08
## AIC: 14136.88
## SBC: 14211.88
## ******************************************************************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")
#15.- CUBIC - P-Splines
# 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")
#16.- REDES NEURONALES
Las redes neuronales (multiperceptrón) también pueden emplearse como funciones smooth (requiere del paquete gamlss.add)
library(gamlss.add)## Warning: package 'gamlss.add' was built under R version 4.2.3## Loading required package: mgcv## This is mgcv 1.8-41. For overview type 'help("mgcv-package")'.##
## Attaching package: 'mgcv'## The following object is masked from 'package:nnet':
##
## multinom## Loading required package: rpartmodelo_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")
Cuando 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.add)
# 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)#17.- INFERENCIA Y PREDICCIÓN
##17.1.- Error en los coeficientes
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##17.2.- 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 sólo aplican si el modelo no incluye términos con funciones smooth. En el caso de que sí los incluya, el resultado obtenido sólo 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)## Single term deletions for
## mu
##
## Model:
## log_recaudacion ~ log_proyecciones
## Df AIC LRT Pr(Chi)
## <none> 13965
## log_proyecciones 1 19904 5940.8 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##17.2.- Intervalo 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
)confint(
modelo_gamlss,
what = c("mu"),
level = 0.95,
robust = TRUE
)## 2.5 % 97.5 %
## mu.(Intercept) 2.9031354 3.3217566
## mu.log_proyecciones 0.8548881 0.8826812confint(
modelo_gamlss,
what = c("sigma"),
level = 0.95,
robust = TRUE
)## 2.5 % 97.5 %
## sigma.(Intercept) 2.0902161 2.2749293
## sigma.log_proyecciones -0.1665666 -0.1508655##17.4.- 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
)predict(
object = modelo_gamlss,
what = "mu",
type = "link",
newdata = data.frame(log_proyecciones = 11)
)## [1] 12.66908predict(
object = modelo_gamlss,
what = "mu",
type = "response",
newdata = data.frame(log_proyecciones = 11)
)## [1] 12.66908predict(
object = modelo_gamlss,
what = "mu",
type = "terms",
newdata = data.frame(log_proyecciones = 11)
)## log_proyecciones
## 1 -0.6808943
## attr(,"constant")
## [1] 13.34997predict(
object = modelo_gamlss,
what = "sigma",
type = "link",
newdata = data.frame(log_proyecciones = 11)
)## [1] 0.4366962predict(
object = modelo_gamlss,
what = "sigma",
type = "response",
newdata = data.frame(log_proyecciones = 11)
)## [1] 1.547586predict(
object = modelo_gamlss,
what = "sigma",
type = "terms",
newdata = data.frame(log_proyecciones = 11)
)## log_proyecciones
## 1 0.1243908
## attr(,"constant")
## [1] 0.3123054predictAll(
object = modelo_gamlss,
type = "response",
newdata = data.frame(log_proyecciones = 11)
)## $mu
## [1] 12.66908
##
## $sigma
## [1] 1.547586
##
## attr(,"family")
## [1] "NO" "Normal"#18.- 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.
K-fold cross validation
El método de K-fold cross validation funciona de la siguiente forma. Las observaciones de entrenamiento se reparten en k conjuntos (folds) del mismo tamaño. El modelo se ajusta con todas las observaciones excepto las del primer fold y se evalúa con las observaciones del fold que ha quedado excluido, obteniendo así la primera métrica. El proceso se repite k veces, excluyendo un fold distinto en cada iteración. Al final del proceso, se generan k valores de la métrica que se agregan (normalmente con la media y la desviación típica), generando la estimación final de validación.
Validación simple
Las observaciones de entrenamiento se reparten en dos conjuntos, uno de entrenamiento y otro de validación/test. El modelo se ajusta con todas las observaciones del conjunto de entrenamiento y se evalúa las observaciones del conjunto de validación, obteniendo así la métrica.
Empleando los datos de rent explicados en el primer de este documento, se entrenan 3 modelos distintos con el objetivo de identificar cuál de ellos predice mejor el precio de alquiler de las viviendas.
modelo_1: predice el precio del alquiler en función de los metros cuadrados de la vivienda, asumiendo que la variable respuesta sigue una distribución normal.
modelo_2: predice el precio del alquiler en función de los metros cuadrados de la vivienda y de su año de construcción, asumiendo que la variable respuesta sigue una distribución normal.
modelo_3: predice el precio del alquiler en función de los metros cuadrados de la vivienda, asumiendo que la variable respuesta sigue una distribución gamma.
modelo_4: predice el precio del alquiler en función de los metros cuadrados de la vivienda y de su año de construcción, asumiendo que la variable respuesta sigue una distribución gamma.
data("rent")
datos <- rent
datos <- datos %>% select(R, Fl, A, H, loc)
# Se nombran 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)## df AIC
## modelo_4 4 28000.20
## modelo_3 3 28107.43
## modelo_2 4 28355.14
## modelo_1 3 28472.35cv_modelo_1 <- gamlssCV(formula = precio ~ metros, family = NO, data = datos,
K.fold = 10, parallel = "multicore", ncpus = 4,
set.seed = 1)## fold 1
## fold 2
## fold 3
## fold 4
## fold 5
## fold 6
## fold 7
## fold 8
## fold 9
## fold 10cv_modelo_2 <- gamlssCV(formula = precio ~ metros + anyo, family = NO, data = datos,
K.fold = 10, parallel = "multicore", ncpus = 4,
set.seed = 1)## fold 1
## fold 2
## fold 3
## fold 4
## fold 5
## fold 6
## fold 7
## fold 8
## fold 9
## fold 10cv_modelo_3 <- gamlssCV(formula = precio ~ metros, family = GA, data = datos,
K.fold = 10, parallel = "multicore", ncpus = 4,
set.seed = 1)## fold 1
## fold 2
## fold 3
## fold 4
## fold 5
## fold 6
## fold 7
## fold 8
## fold 9
## fold 10cv_modelo_4 <- gamlssCV(formula = precio ~ metros + anyo, family = GA, data = datos,
K.fold = 10, parallel = "multicore", ncpus = 4,
set.seed = 1)## fold 1
## fold 2
## fold 3
## fold 4
## fold 5
## fold 6
## fold 7
## fold 8
## fold 9
## fold 10CV(cv_modelo_1, cv_modelo_2, cv_modelo_3, cv_modelo_4)## val[o.val]
## cv_modelo_4 27998.08
## cv_modelo_3 28105.69
## cv_modelo_2 28353.58
## cv_modelo_1 28470.61# 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)## Pred.GD
## test_modelo_4 5569.789
## test_modelo_3 5598.020
## test_modelo_2 5623.987
## test_modelo_1 5656.076modelo_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 28019
## + 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 27809
## <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
## - anyo 1 27818
## ---------------------------------------------------
## 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.8modelo_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 28577
## + 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 28401mejor_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##INFO DE LA SESIÓN
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)## # A tibble: 101 × 3
## package loadedversion source
## <chr> <chr> <chr>
## 1 abind 1.4-5 CRAN (R 4.2.0)
## 2 backports 1.4.1 CRAN (R 4.2.0)
## 3 base64enc 0.1-3 CRAN (R 4.2.0)
## 4 broom 1.0.5 CRAN (R 4.2.3)
## 5 bslib 0.4.2 CRAN (R 4.2.2)
## 6 cachem 1.0.6 CRAN (R 4.2.2)
## 7 callr 3.7.3 CRAN (R 4.2.3)
## 8 car 3.1-2 CRAN (R 4.2.3)
## 9 carData 3.0-5 CRAN (R 4.2.3)
## 10 cli 3.6.0 CRAN (R 4.2.2)
## # ℹ 91 more rows