Implementando distribuciones en gamlss
Introducción
GAMLSS (Generalized Additive Models for Location, Scale and Shape) son una poderosa evolución de los modelos estadísticos clásicos, porque no se limitan a explicar solo la media de una variable, sino que permiten modelar también su variabilidad, asimetría y forma completa de la distribución. En lugar de encasillar los datos en distribuciones rígidas, GAMLSS abre la puerta a más de 100 distribuciones flexibles, ofreciendo un marco donde los datos se ajustan a la realidad y no al revés. Esto convierte a GAMLSS en una herramienta moderna, versátil y apasionante para quienes buscan ir más allá de los modelos tradicionales y explorar el comportamiento completo de los datos.
En esta publicación se muestran los pasos necesarios para implementar una distribución en gamlss y aprovechar toda la infraestructura de los paquetes gamlss y gamlss2.
En la siguiente lista de reproducción se muestran videos cortos con una explicación de lo que son los modelos gamlss.
Pasos
1. Identificar la distribución.
Lo primero que se debe hacer es buscar en un artículo científico o libro una distribución estadística interesante (discreta o continua) que no esté implementada en gamlss. Para esa distribución debemos elegir un nombre corto que identifique la distribución, algunos nombres cortos para distribuciones estadísticas que ya están en gamlss son: NO para la normal, BE para la beta y PO para la Poisson.
En los siguientes enlaces hay un listado de las distribuciones que ya están implementadas en gamlss:
A partir de aquí vamos a suponer que el nombre corto que elegimos para la distribución que queremos implementar es XXX.
2. Identificar el número de parámetros.
Lo siguiente es identificar el número de parámetros de la distribución XXX seleccionada. Si la distribución tiene 4 o menos parámetros, la distribución se puede implementar en gamlss.
3. Identificar el rol de los parámetros.
Debemos identificar el rol de cada uno de los parámetros en la distribución XXX, es decir, identificar cuál parámetro está asociado con la localización (media, mediana), cuál con variabilidad (varianza o desviación), cuál con la asimetría y cuál con la kurtosis.
4. Cambiar los nombres de los parámetros.
Las distribuciones en gamlss deben tener máximo 4 parámetros y se deben llamar \mu, \sigma, \nu y \tau, no es posible usar otras letras griegas para nombrar a los parámetros.
La convención para hacer el cambio de nombres es la siguiente:
- El parámetro asociado con la localización llamarlo \mu,
- El parámetro asociado con la variabilidad llamarlo \sigma,
- Los parámetros de asimetría y kurtosis llamarlos \nu y \tau.
Si la localización, varianza, asimetría y kurtosis de la distribución XXX no están relacionadas con un solo parámetro, podemos hacer el cambio de nombres de cualquier forma.
Si la distribución seleccionada tiene un sólo parámetro, use \mu; si tiene dos, use \mu y \sigma; si tiene tres, use \mu, \sigma y \nu; si tiene cuatro, use \mu, \sigma, \nu y \tau.
5. Construcción de las funciones dXXX, pXXX, qXXX y rXXX.
El siguiente paso es construir las funciones dXXX, pXXX, qXXX y rXXX para la distribución XXX de interés. A continuación se muestran varios videos en los cuales se explica cómo construir estas funciones para varios casos.
📌 : en los videos y explicaciones que siguen, vamos a usar como ejemplo una distribución con dos parámetros (\mu y \sigma) por comodidad, sin embargo, las explicaciones se extiden de manera similar a distribuciones con 3 y 4 parámetros.
Caso: distribución discreta
Función de masa dXXX.
Función probabilidad acumulada pXXX.
Función de cuantiles qXXX.
Función generadora de números aleatorios rXXX.
Casi todas las distribuciones propuestas en la literatura cuentan con una expresión matemática para su función de masa de probabilidad f(x), sin embargo, algunas distribuciones no cuentan una expresión explícita para la función acumulada F(x) o para la función de cuantiles Q(p). En estos casos es necesario recurrir a la definición de F(X) y de Q(p) para poder crear las funciones pXXX y qXXX. En los siguientes dos videos se muestra cómo crear estas funciones.
Función probabilidad acumulada pXXX cuando no se tiene fórmula para F(X).
Función de cuantiles qXXX cuando no se tiene fórmula para Q(p).
Caso: distribución continua
⏰ Pronto vamos a tener aquí otros videos sobre cómo crear las siguientes funciones para el caso de una distribución continua.
Función de densidad dXXX.
Función probabilidad acumulada pXXX.
Función de cuantiles qXXX.
Función generadora de números aleatorios rXXX.
6. Revisando la función dXXX(x)
Se debe revisar siempre que las funciones que se han creado estén haciendo lo correcto. A continuación se explica una prueba sencilla para revisar la función dXXX(x) .
- Si la distribución es discreta: calcule la sumatoria de las probabilidades obtenidas con
dXXX(x)para un vectorxcon muchos valores posibles de la distribución. A continuación un ejemplo para un vectorxcon mil valores consecutivos y dos valores arbitrários de \mu y \sigma:
x <- 0:999
sum(dXXX(x=x, mu=1.7, sigma=2.3)) Lo que tiene que obtener del código anterior es que la suma sea muy cercana a 1.
- Si la distribución es continua: calcule la integral de
dXXX(x)desde un valor muy a la izquierda hasta un valor muy a la derecha, por ejemplo, desde -99 hasta 99. A continuación un ejemplo de la integral para dos valores arbitrários de \mu y \sigma:
integrate(dXXX, lower=-99, upper=99, mu=1.7, sigma=2.3) Lo que debe obtener es que la integral sea muy cercana a 1.
7. Derivadas de \log(f(x))
El siguiente paso consiste en escribir la función \log(f(x)) y calcular las primeras derivadas parciales con respecto a cada uno de los parámetros de la distribución.
Para explicar este paso vamos a utilizar como ejemplo la distribución DPERKS discreta de dos parámetros que fue explicada en este enlace.
A continuación la función f(x) para la distribución DPERKS.
Lo siguiente es calcular \log(f(x)) y calcular las primeras derivadas de \log(f(x)) con respecto a los dos parámetros para obtener dldm y dldd así:
8. Crear funciones de R para las primeras derivadas
Luego de haber calculado analíticamente las primeras derivadas es necesario programarlas en R.
Como ejemplo vamos a programar dldm para la distribución DPERKS.
dldm_manual = function(y, mu, sigma) {
part1 <- 1/mu + 1/(mu+1)
part2 <- - exp(sigma*y)/(1+mu*exp(sigma*y))
part2 <- - exp(sigma*(y+1))/(1+mu*exp(sigma*(y+1)))
dldm <- part1 + part2 + part3
dldm
}📌 : al programar dldm debemos usar la letra y en lugar de la x.
✳️ : se debe crear una función por cada primer derivada. Si la distribución tiene tres parámetros (\mu, \sigma, \nu), se deben crear tres funciones.
9. Verificar que las derivadas estén correctas
Luego de crear cada derivada manual es necesario verificar que esté correcta.
En el siguiente video se muestra cómo hacer la verificación utilizando como ejemplo la distribución GEG.
10. Segundas derivadas de \log(f(x))
Un elemento necesario para crear la familia es d2ldm2, d2ldd2 y d2ldmdd. Estos elementos se pueden obtener de dos formas alternativas:
- como el valor esperado de la segunda derivada parcial,
- como menos el producto de la primera derivada parcial.
Siempre se debe intentar la primera forma, en caso de no llegar a una expresión a la cual se le pueda obtener el valor esperado, se puede usar la segunda forma.
A continuación se muestran dos ejemplos, en ninguno de ellos se llega a una expresión sencilla a la cual se le aplica el valor esperado E(\cdot).
Necesitamos obtener d2ldd2 así:
Necesitamos obtener d2ldmdd así:
📌 : para los elementos d2ldm2, d2ldd2 y d2ldmdd no es necesario hacer chequeos computacionales de las derivadas.
11. Creando la función familia
12. Verificando la familia
Este paso es muy importante para saber si la función de familia quedó bien construída.
Supongamos que la distribución de interés es la DPERKS para la cual ya tenemos construidas las funciones dDPERKS, pDPERKS, qDPERKS, rDPERKS y DPERKS.
El código que se muestra a continuación sirve para crear una muestra aleatoria de n valores de nuestra distribución con valores apropiados de \mu y \sigma, luego con la ayuda de la función gamlss ajustamos un modelo y luego solicitamos los parámetros estimados.
n <- 5000
# True parameters are:
true_mu <- 2.5
true_si <- 4.4
y <- rDPERKS(n=n, mu=true_mu, sigma=true_si)
library(gamlss)
mod <- gamlss(y ~ 1, family=DPERKS)
coef(mod, what="mu")
coef(mod, what="sigma")Taller práctico
Para aplicar lo aprendido en los videos creamos dos talleres prácticos para que usted implemente una distribución en gamlss.
