En el estudio y modelamiento de datos en diversas áreas del saber los Modelos Lineales Generalizados, en adelante GLM, son unas de las técnicas más utilizadas, difundidas y estudiadas. Los GLM fueron introducidos por Nelder and Wedderburn (1972) y, de manera general, pueden ser entendidos como una extensión de las técnicas de regresión lineal a la familia exponencial uniparamétrica, de ahí el nombre generalizados (Gelfand and Dalal (1990); Dey, Gelfand, and Peng (1997)).
En su planteamiento, los GLM cuentan con tres componentes: primero, una componente aleatoria que es la variable a modelar. Segundo, una componente sistemática, que describe los cambios en la variable aleatoria como una función de las covariables. Finalmente, la tercera componente es la función de enlace, la cual conecta la componente aleatoria con la componente sistemática del modelo (McCullagh and Nelder (2019)).
Debido a su amplio uso, existe una extensión de los GLM muy útil en contextos de dispersión variable ya que le asigna un predictor lineal al parámetro de dispersión. A esta técnica se le ha dado el nombre de Modelos Lineales Doblemente Generalizados (DGLM). Los DGLM pueden ser entendidos como una generalización de los GLM a la familia exponencial biparamétrica en donde se le asigna un predictor lineal a cada uno de los parámetros. La diferencia entre los GLM y los DGLM es que los segundos buscan modelar no sólo la relación de un parámetro, usualmente la media, con un conjunto de covariables, sino que buscan modelar la relación de dos parámetros, usualmente la media y el parámetro de dispersión con un conjunto de covariables (Smyth and Verbyla (1999)).
En Smyth and Verbyla (1999) se plantea que, en el caso en el que los conjuntos de datos presenten dispersión variable, una manera de afrontarla consiste en llevar a cabo el modelamiento de esta; es decir, aplicar una estructura de regresión no solo sobre la media sino también sobre la dispersión del conjunto de datos.
Los DGML pueden ser entendidos como una generalización de los GLM para el caso en el que existe dispersión variable, por lo tanto, acorde con [Smyth & Verbyla, 1999] puede ser entendido como un conjunto de dos submodelos en el cual por un lado se modela la media de la variable y por el otro se modela la dispersión o algún otro parámetro de interés.
Sean \(Y_i, i=1,2,...,n\) variables aleatorias independientes. Los modelos lineales doblemente generalizados, DGLM, están definidos por las siguientes tres componentes:
La primera componente del modelo, la componente aleatoria, es el vector \(Y = (y_1,..., y_n)\) con elementos independientes con distribución en la familia exponencial biparamétrica. Por independencia se tiene que \(Cov(y_i, y_j) = 0\); sin embargo, como los elementos no son idénticamente distribuidos se puede atribuir a cada uno de ellos \(E(y_i) = \mu_i\) y \(Var(y_i) = \sigma^2_i\).
La segunda componente del modelo son las componentes sistemáticas, conformadas por los predictores lineales que componen el vector \(\eta = (\eta_1, \eta_2)\) y se define \(\eta_1 = X\beta\) y \(\eta_2 = Z\gamma\) donde X es una matriz de dimensión \(n \times p\) y Z es una matriz de dimensión \(n\times l\), \(\beta\) es un vector \(p \times 1\) y \(\gamma\) es un vector \(l \times 1\). Las matrices X y Z contienen en cada una de sus filas los valores de las covariables de cada uno de las n observaciones. Las variables en X modelarán un componente, usualmente la media; y las variables Z modelarán otra componente, usualmente la dispersión.
La tercera componente del modelo son las funciones de enlace, las cuales, en este caso, conectan las componentes sistemáticas con las componentes aleatorias, es decir, \(h_1(\theta) = \eta_1\) y \(h_2(\tau) = \eta_2\). En este caso se está modelando el comportamiento de los parámetros \(\theta\) y \(\tau\), además, \(h_1\) y \(h_2\) son funciones monótonas y diferenciables.
Descripción: el paquete dglm (Smyth (1989)) sirve para ajustar un modelo lineal generalizado para la media así como también un modelo adicional de enlace para la dispersión.
Devuelve un objeto de clase “dglm”, que hereda de “glm” y “lm”.
Un objeto “dglm” consiste de un objeto “glm” con los siguientes componentes adicionales:
´dispersion.fitthe´ submodelo de la dispersión: un objeto “glm” que representa el modelo ajustado para las dispersiones. Las respuestas para este modelo son los componentes de desviación del modelo lineal generalizado original. Los pesos anteriores son 1 y la dispersión o escala de este modelo es 2.
El componente ´iterthis´ ahora representa el número de iteraciones externas utilizadas para ajustar los modelos acoplados de dispersión y de media. En cada iteración externa, se realiza una IRLS para cada uno de los submodelos de media y dispersión.
´methodfitting´ método usado: “ML” si se utilizó máxima verosimilitud o “REML” si se utilizó máxima verosimilitud restringida.
´m2loklikminus´ el doble de la log-verosimilitud o verosimilitud perfilada del modelo ajustado.
Carga de los paquetes necesarios para la ejecución del modelo lineal doblemente generalizado:
library(dglm)
## Loading required package: statmod
library(faraway)
A continuación, trabajaremos con la base de datos motorins (Hallin and Ingenbleek (1983)) del paquete faraway la cual consiste de un data-frame con 1797 observaciones de reclamaciones de seguros de automóviles de terceros en Suecia en 1977 de las siguientes 8 variables :
Con las siguientes instrucciones, invocamos la BD motorins y transformamos a factor las variables necesarias:
attach(motorins)
head(motorins)
## Kilometres Zone Bonus Make Insured Claims Payment perd
## 1 1 1 1 1 455.13 108 392491 3634.176
## 2 1 1 1 2 69.17 19 46221 2432.684
## 3 1 1 1 3 72.88 13 15694 1207.231
## 4 1 1 1 4 1292.39 124 422201 3404.847
## 5 1 1 1 5 191.01 40 119373 2984.325
## 6 1 1 1 6 477.66 57 170913 2998.474
Kilom <- as.factor(Kilometres)
Bonus <- as.factor(Bonus)
Make <- as.factor(Make)
Ajustamos un modelo Poisson tipo “dglm” donde la variable respuesta es Claims, además de modelar la dispersión con una sola covariable, Make, como sigue:
dglm.claims <- dglm(formula = Claims~Kilom+Bonus+Make,
offset=log(motorins$Insured), ~Make,
family=poisson(link="log"), data=motorins)
Y vemos su resumen:
summary(dglm.claims)
##
## Call: dglm(formula = Claims ~ Kilom + Bonus + Make, dformula = ~Make,
## family = poisson(link = "log"), data = motorins, offset = log(motorins$Insured))
##
## Mean Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.53686410 0.09119122 27.819169 2.543592e-170
## Kilom.L -0.89636078 0.06342735 -14.132085 2.409036e-45
## Kilom.Q -0.33345159 0.05923268 -5.629520 1.807115e-08
## Kilom.C 0.40989671 0.05966601 6.869853 6.426817e-12
## Kilom^4 0.09150786 0.05330771 1.716597 8.605282e-02
## Bonus 0.27818635 0.01268013 21.938758 1.108873e-106
## Make2 -1.37718576 0.09046589 -15.223261 2.478636e-52
## Make3 -1.63895661 0.09688573 -16.916387 3.406984e-64
## Make4 -1.48381719 0.13063045 -11.358892 6.698885e-30
## Make5 -1.23695463 0.08887854 -13.917359 4.969268e-44
## Make6 -0.85680194 0.09250011 -9.262713 1.993024e-20
## Make7 -1.56775096 0.08931249 -17.553547 5.587902e-69
## Make8 -2.15628348 0.09820706 -21.956502 7.506011e-107
## Make9 1.93482328 0.09062734 21.349223 3.964744e-101
## (Dispersion Parameters for poisson family estimated as below )
##
## Scaled Null Deviance: 4192.811 on 1796 degrees of freedom
## Scaled Residual Deviance: 1796.999 on 1783 degrees of freedom
##
## Dispersion Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.8280955 0.09386465 40.783143 0.000000e+00
## Make2 -1.3763424 0.13625957 -10.100886 5.474538e-24
## Make3 -1.5244284 0.14158576 -10.766820 4.937633e-27
## Make4 -0.5337588 0.14735603 -3.622239 2.920637e-04
## Make5 -1.3268728 0.13608567 -9.750276 1.839680e-22
## Make6 -0.7623476 0.13507233 -5.643996 1.661484e-08
## Make7 -1.6575479 0.13770569 -12.036887 2.273742e-33
## Make8 -1.9927300 0.14272795 -13.961737 2.668514e-44
## Make9 2.0493876 0.13041202 15.714714 1.199209e-55
## (Dispersion parameter for Gamma family taken to be 2 )
##
## Scaled Null Deviance: 4490.145 on 1796 degrees of freedom
## Scaled Residual Deviance: 2642.288 on 1788 degrees of freedom
##
## Minus Twice the Log-Likelihood: 14997.46
## Number of Alternating Iterations: 6
Usando los resultados del resumen anterior, podemos escribir el modelo como sigue:
\[Claims \thicksim Poisson (\widehat\mu_i),\] \[log(\widehat\mu_i) = 2.54 - 0.90Kilom.L + \dotsb + 1.94Make9,\] \[log(\widehat\phi_i)=3.83+1.38Make2+ \dotsb + 2.05Make9\]
Y, por último, vemos gráficos que nos ayudan a validar el ajuste del modelo a partir de los residuales:
par(mfrow=c(2,2))
plot(dglm.claims)
Se logra evidenciar un comportamiento aceptable en los gráficos de residuales aún cuando algunos de estos se alejan de la línea de 45° en el Q-Q Plot. Las anteriores gráficos muestran un buen diagnóstico del modelo ajustado.
Nota: Para conocer otras publicaciones relacionadas con glm visite este enlace.