1 Algoritmos en gamlss

Para la obtenciĆ³n de los parĆ”metros en modelos gamlss estĆ”n disponibles tres algoritmos:

  • RS(): The default method is the RS algorithm, which does not requires accurate starting values for Āµ, Ļƒ, Ī½ and Ļ„ to ensure convergence (the default starting values, often constants, are usually adequate). This method is faster for larger data sets.
  • CG(): The CG algorithm, which can be better for distributions with potentially highly correlated parameter estimates.
  • mixed(): This is a mixture of the above two algorithms which starts with the RS algorithm and finishes with the CG.

2 ComparaciĆ³n

Vamos a ajustar el mismo modelo para explicar la media \(\mu\) y la desviaciĆ³n \(\sigma\) de la variable y en funciĆ³n de x usando un modelo cuadrĆ”tico. Lo primero es cargar el paquete y disponer de los datos.

library(gamlss)
data(abdom)
with(abdom, plot(x=x, y=y, las=1))

A continuaciĆ³n vamos a ajustar tres modelos que solo difieren en el algoritmo o mĆ©todo.

m1 <- gamlss(y ~ x + I(x^2), sigma.fo= ~x + I(x^2), 
             family=NO, data=abdom)
## GAMLSS-RS iteration 1: Global Deviance = 4793.639 
## GAMLSS-RS iteration 2: Global Deviance = 4792.8 
## GAMLSS-RS iteration 3: Global Deviance = 4792.8
m2 <- gamlss(y ~ x + I(x^2), sigma.fo= ~x + I(x^2), 
             family=NO, data=abdom, method=CG())
## GAMLSS-CG iteration 1: Global Deviance = 6166.027 
## GAMLSS-CG iteration 2: Global Deviance = 5635.653 
## GAMLSS-CG iteration 3: Global Deviance = 5208.609 
## GAMLSS-CG iteration 4: Global Deviance = 4942.862 
## GAMLSS-CG iteration 5: Global Deviance = 4827.348 
## GAMLSS-CG iteration 6: Global Deviance = 4795.882 
## GAMLSS-CG iteration 7: Global Deviance = 4792.856 
## GAMLSS-CG iteration 8: Global Deviance = 4792.801 
## GAMLSS-CG iteration 9: Global Deviance = 4792.8 
## GAMLSS-CG iteration 10: Global Deviance = 4792.8
m3 <- gamlss(y ~ x + I(x^2), sigma.fo= ~x + I(x^2), 
             family=NO, data=abdom, method=mixed(2, 20))
## GAMLSS-RS iteration 1: Global Deviance = 4793.639 
## GAMLSS-RS iteration 2: Global Deviance = 4792.8 
## GAMLSS-CG iteration 1: Global Deviance = 4792.8

De los resultados anteriores vemos que m1 converge en 3 iteraciones, que m2 demora hasta la iteraciĆ³n 10 para la convergencia y que m3 convege tambiĆ©n rĆ”pido porque inicia con 2 iteraciones de RS y luego se mejoran las estimaciones con una iteraciĆ³n de CG.

Para extraer y comparar los coeficientes de los modelos vamos a crear la siguiente funciĆ³n.

my_coefs <- function(mod, digits=6) {
  res <- c(coef(mod, what="mu"), coef(mod, what="sigma"))
  return(round(res, digits=digits))
}

Ahora hacemos la comparaciĆ³n asĆ­:

cbind(my_coefs(m1), my_coefs(m2), my_coefs(m3))
##                   [,1]       [,2]       [,3]
## (Intercept) -96.922008 -96.915642 -96.922008
## x            13.672570  13.671991  13.672570
## I(x^2)       -0.060315  -0.060303  -0.060315
## (Intercept)   1.403001   1.403485   1.402300
## x             0.038545   0.038506   0.038603
## I(x^2)        0.000068   0.000069   0.000067

De la salida anterior vemos que los coeficientes estimados estƔn muy cercanos.