UNIVERSIDAD NACIONAL DE LA PLATA

FACULTAD DE CIENCIAS AGRARIAS Y

FORESTALES

CURSO: Modelos Lineales Generalizados Aplicados

a las Ciencias Biológicas

MATERIAL DE APOYO
Notas sobre Modelos Lineales Generalizados:
Una Introducción

Prof.: Dra. María del Pilar Díaz (Universidad Nacional de Córdoba).

PRESENTACIÓN:

Durante muchos años el área de modelación estuvo restringida a la de los

modelos normales clásicos, los cuales eran utilizados para describir la
mayoría de los fenómenos aleatorios, aún en los casos en que era bastante
razonable suponer distribución no normal para el comportamiento de la
variable en estudio. En algunos de estos casos, y aún en la actualidad, se
utilizan transformaciones para intentar lograr normalidad. Las transfor-
maciones deben garantizar también simultáneamente los otros supuestos
asumidos en un modelo clásico, como constancia de la varianza y linealidad
(aditividad) en los parámetros, lo cual ocurre raramente. Así, el uso de
datos transformados como base del análisis estadístico es adecuado sólo
cuando la escala que hace que los supuestos requeridos se cumplan, tiene
significado en el área de estudio, dado que las conclusioes se aplican a las
poblaciones transformadas (Mead et al, 1993). Nótese que mientras en
supuesto de aditividad concierne a la construcción del modelo, los otros
dos, normalidad y homogeneidad de varianzas se relacionan con la variación
aleatoria y no hay garantía de que ambos requisitos resulten de la misma
transformación (Mead et al, 1993; Box y Cox, 1982).

Una propuesta interesante, basada en un concepto innovador de

modelación, surgió en 1972, con el artículo de Nelder y Wedderburn sobre

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
los modelos lineales generalizados (MLG). La idea básica consiste en abrir
el abanico de posibilidades u opciones para la distribución de la variable
respuesta (esto es, “relajar” el supuesto de distribución normal), siempre
y cuando pertenezca a una familia más amplia de distribuciones: la familia
exponencial, así como permitir que la relación entre el valor esperado
(media) de la variable y la combinación lineal de los parámetros (parte
sistemática del modelo) no sea siempre la identidad, sino cualquier función
monótona.

Por ejemplo, para datos de conteo, en vez de aplicar la transformación

y para lograr normalidad, podemos suponer que la distribución para la
variable es Poisson (ya que es una variable discreta) y que la relación entre
el valor esperado y la combinación lineal en los parámetros es
log µ = η = β 0 + β 1 x1 + β 2 x 2 + L + β p −1 x p −1 , la cual asegura que cualesquiera
sean los valores de los parámetros siempre dará un valor positivo para µ
(el conteo esperado). Podemos proceder de manera semejante en el caso
de los datos sean proporciones y pensar que la distribución Binomial es una
candidata razonable para describir el comportamiento de la variable
aleatoria proporción de éxitos. Una relación funcional entre µ, la
proporción esperada de éxitos, y η, el predictor lineal, puede ser del tipo
µ
log , llamada función logística.
1− µ

Los MLG permiten trabajar con variables de tipo continuas, discretas y

categóricas, pertenecientes a la familia exponencial (funciones de
densidad normal, gama, normal inversa, distribuciones de probabilidad
Poisson, Binomial, Multinomial, entre otras). Por consiguiente, el supuesto
de varianzas constantes puede ser relajado, ya que se permite usar
funciones de varianza, por ejemplo, en una Binomial, la varianza nµ(1-µ) es
función de la media µ. Por otra parte la aditividad de los efectos
sistemáticos puede ocurrir en la escala de una transformación monótona
de la media. Es entonces clara la ventaja de los MLG frente a la
transformación de datos.

Como muchas de estas relaciones funcionales pueden ser del tipo no

lineal, Nelder y Wedderburn (1972) propusieron la estimación de los
parámetros a través de un proceso iterativo, a la vez que presentaron una
medida de bondad de ajuste (ya no una suma de cuadrados, como es en el
caso de un modelo lineal normal), llamanda deviance, que es utilizada en el

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
ajuste del modelo y en las etapas de diagnóstico. Como en rigor el método
de estimación se basa en el principio de máxima verosimilitud, los
estimadores que se obtienen tienen buenas propiedades estadísticas.

A partir de trabajo de Nelder y Wedderburn, en 1972, fueron publicados

innumerables trabajos científicos sobre MLG. Mas, fue a partir del
desarrollo del primer software para la estimación de los mismos que este
marco metodológico cobró importancia, principalmente en el campo de las
Ciencias Biológicas. Ese sistema estadístico, llamado GLIM (Generalized
Linear Interactive Models1, (vide Aitkin et al., 1989, Francis et al. 1994),
fue desarrollado exclusivamente para el ajuste de MLG, y fue por mucho
tiempo, única alterativa para el análisis de datos. Hoy existen otros
sistemas, como el S-Plus2, R3, SAS4, STATA5, por citar algunos.

Esta teoría ha sido vastamente extendida ya. Wedderburn (1974)

proporcionó la base teórica para los modelos de cuasi-verosimilitud, los
cuales generalizan los MLG a situaciones más generales, incluyendo datos
correlacionados. Jørgensen (1983), con la construcción de los modelos de
dispersión, amplia aún más uno de los aspectos centrales de la definición
de los MLG: la distribución de la variable respuesta.

En 1986, Liang y Zeger (1986) extienden los modelos de cuasi-

verosmilitud presentando las ecuaciones generalizadas de estimación
(conocidas como GEE), las cuales permiten el abordaje de estudio
longitudinales (tiempo o espacio) descriptos por variables aleatorias no
normales correlacionadas.

Por otro lado, Hastie y Tibshirani (1990) presentam los modelos aditivos
generalizados (GAM), que suponen un predictor lineal que puede ser
formado por funciones semiparamétricas, adecuadas para descripciones
de patrones no lineales que requieren de suavizados. Breslow y Clayton
(1993) fueron los primeros en constriur el marco teórico para los modelos
lineales generalizados mixtos (GLMM), en el sentido de admitir la inclusión
de efectos aleatórios (normales) en el predictor lineal. Muchos de esos
resultados se discuten en McCulloch y Searle (2001). Actualmente las

1
(http://www.nag.co.uk/stats/gdge_soft.asp),
2
(http://www.insightful.com)
3
(http://www.r-project.org)
4
(http://www.sas.com),
5
(http://www.stata.com),

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
aplicaciones de los MLG pueden encontrarse en casi todas las disciplinas
científicas, siendo un libro óptimo de referencia el McCullagh y Nelder
(1989), que será abordado en este curso.

ALGUNAS IDEAS...1.

1. Modelo lineal generalizado

Sean Y1 ,..., Yn variables aleatorias (de respuestas) independientes, que
originarán nuestro conjunto de observaciones. Un MLG es un modelo que
vincula características de estas respuestas (variables “dependientes”) con
otras variables “explicativas”. Tenemos que considerar tres componentes:

1. La componente aleatoria (la distribución de las Yi ). En general, se

supone que las Yi son independientes, con una distribución que sea
una familia exponencial lineal (uniparamétrica).
2. La componente sistemática o predictor lineal, que indica la relación
entre las covariables (variables explicativas). Este es un modelo
lineal (es decir, los parámetros entran linealmente al modelo). Por
ej., α + β1 x1 + β 2 x2 .
3. La función de enlace o ligamiento, que es la que vincula los valores
esperados de las Yi con la componente sistemática. Por ejemplo,
g ( µi ) = log( µi ) = α + β1 x1i + β 2 x2i .

Es sólo a través de la función de enlace aplicada al predictor lineal que las

covariables ejercen efecto sobre la variable respuesta. El enlace canónico
(una posibilidad que se deduce desde la formulación dentro de la familia
exponencial) se torna una escala adecuada en la cual el valor esperado se
puede describir mediante una función lineal en los parámetros. McCullagh
y Nelder (1989) plantean que la elección de la escala es un aspecto
importante en la selección del modelo. La escala adecuada para el análisis
clásico de regresión lineal, por ejemplo, es aquella que combina varianza
constante, errores independientes distribuidos aproximadamente
normales y aditividad en los efectos sistemáticos. Sin embargo, con la
introducción de los MLG la normalidad y la varianza constante ya no son
necesarias y la aditividad de los efectos, aunque importante, puede
lograrse en otra escala. Lo que debe conocerse es la manera en que la
varianza depende de la media, esto es, la relación entre la media y la
varianza establecida por el modelo probabilístico seleccionado, mientras

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
que la aditividad se postula como propiedad de una cierta función de la
esperanza de la variable respuesta.
Algunos ejemplos de modelos lineales generalizados pueden ser:
a. Yi  (α + β1 x1i + β 2 x2i , σ 2 ) ;
recordando que denotamos Y N(mu, signma2), observamos que el valor
esperado es exactamente una función lineal en los parámetros;
b. Yi Poisson( µi ); log( µi ) = α + τ i ,
obsérvese que el valor esperado no es, como sucede en el caso a), una
función monótona de una función lineal en los parámetros.
 π 
c. Yi Bin(π i ,  ); log  i  = β 0 + β1 xi ,
 1−πi 
similar al caso b), una función no lineal del valor esperado es igual a una
función lineal de los parámetros.

2. Modelos más comunes:

2.1 Regresión logística. Interpretación.

Para modelar datos binarios (dos categorías) que provengan de variables

aleatorias distribuídas según una distribución binomial, podemos
considerar modelos que relacionen el parámetro π (que es la proporción
esperada de “éxitos”) con la(s) variable(s) explicativas (mal llamadas
“independientes”). El predictor lineal más simple sería π i = α + β xi , el cual
define a la función de enlace como la identidad. Este modelo que iguala la
probabilidad de éxito a una función lineal (llamado algunas veces modelo de
probabilidad lineal) tiene algunos problemas serios en su aplicabilidad: los
valores predichos no están acotados, dado que el término de la derecha
puede tomar cualquier valor en los ℜ, aunque es óbvio que deberían estar
siempre entre 0 y 1.

El modelo de regresión logística usa como función de enlace el “logit”, lo

que hace que los valores predichos queden siempre entre 0 y 1:
 π ( xi )  exp(α + β xi )
log   = α + β xi ; π ( xi ) = .
 1 − π ( xi )  1 + exp(α + β xi )
Observar que hemos puesto como componente sistemática una regresión
lineal simple, pero podríamos haber puesto cualquier otro modelo lineal.
Por ejemplo, supongamos que los parámetros sean α = −4, β = 0.5 . Luego, se
tiene el comportamiento para el valor esperado como describe Fig. 1.

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
¿Cómo interpretamos el coeficiente de regresión β ?
• tasa de incremento: a medida que | β | crece, π ( xi ) cambia más
rápidamente a medida que x se incrementa,

Ecuación logística

1.0

0.9

0.8

0.7

0.6
pi(x)

0.5

0.4

0.3

0.2

0.1

0.0
0 5 10 15 20

Fig.1
• la pendiente de la curva en cualquier punto es βπ (1 − π ) , y la
pendiente es máxima cuando π = 0.5 .
• el valor de x cuando la pendiente es máxima (también el punto de
inflexión) es x = − α
β (también se llama valor efectivo mediano o
dosis letal 50%, denotado comunmente por LD50),
• el valor exp( β ) representa el cociente de chances (odds-ratio)
cuando las x aumentan en una unidad.

La interpretación como cociente de chances es la que hace muy atractiva

esta función de enlace (que se prefiere a otras como probit o
complemento log-log). La propiedad simétrica del cociente de chances nos
permite usarlo tanto en estudios prospectivos como en retrospectivos,
(caso-control, etc.), y por lo se torma muy útil en estudios de
Epidemiología, Sociología, Ecología.

La estimación de los parámetros del modelo se realiza por máxima

verosimilitud, de manera que para la inferencia podemos usar todas las
propiedades de los estimadores máximo-verosímiles. En general, las
propiedades son asintóticas, por lo que para que se cumplan,
aproximadamente, bien en una situación particular necesitaremos tamaños
de muestra grandes. Computacionalmente, en los MLG, los estimadores se

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
obtienen usando el algoritmo de Newton-Raphson (o algún refinamiento de
este algoritmo).

2.2 Inferencia

Para verificar hipótesis existen tres métodos basados en la función de

verosimilitud: la prueba de Wald, la prueba del cociente de verosimilitud y
la prueba score (ya conocida directamente como prueba escore).
Asintóticamente las tres pruebas son equivalentes, pero para muestras no
demasiado grandes los resultados pueden ser bastante diferentes (p.ej.,
la prueba de Wald tiende a ser más liberal).

Para ver la diferencia entre los tres métodos, vamos a considerar el caso
más simple posible: probar β = 0 en un GLM que tiene ese único parámetro
(es decir, la función de verosimilitud es función sólo de β ). Para estimar
el parámetro, maximizamos la función de verosimilitud, o su logaritmo,
L( β ) . El estimador será denotado como β̂ , y el valor máximo de la log-
verosimilitud es L( βˆ ) . Veamos cómo se formulan las tres pruebas:
2
 βˆ 
• La prueba de Wald es basada en el estadístico χ = 2
, es decir
Wald
s.e.βˆ 
 
el valor estimado dividido por su error estándar (asintótico). Este
error estándar asintótico se calcula a partir de la curvatura de la
función de log-verosimilitud en su máximo. Así, intuitivamente, si la
curva es muy “amplia” o “abierta”, el error estándar es grande (la
estimación no es muy precisa), si la curva es muy “cerrada”, el error
estándar es pequeño (la estimación es precisa). Esta prueba trata de
comprobar si la diferencia (en el eje horizontal) entre el valor
estimado y el valor hipotetizado en la nula es suficientemente pequeña.

• La prueba de verosimilitud usa como estadístico a χ LR

2
(
= −2 L(0) − L( βˆ ) , )
donde los L denotan a la función de verosimilitud evaluada en su valor
medio (valor nulo) y en el estimador que se propone como de máxima
verosimilitud β̂ . Esta prueba observa si la diferencia (en el eje
vertical) entre el máximo obtenible de la log-verosimilitud y el valor de
la log-verosimilitud cuando la hipótesis nula es cierta es
suficientemente pequeña.

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
 2
 L '(0) 
• El estadístico para la prueba del escore es χ = 2
Score  . Esto es,
 s.e.L '(0) 
evalúa si la derivada en el valor postulado es lo suficientemente
cercana a 0. Debemos recordar que para la mayoría de los modelos
estudiados, la función de verosimilitud es cóncava (al menos cerca del
máximo), y que por lo tanto a medida que la derivada se acerca más a 0,
significa que la log-verosimilitud, bajo la hipótesis nula, está más cerca
del máximo.

La prueba de Wald es la más simple computacionalmente, pero es la que

menos se recomienda, ya que puede ser excesivamente liberal cuando los
tamaños de muestra son pequeños y no es invariante ante cambios de
parametrización. Las otras dos pruebas, basados en los estadísticos
escore y cociente de verosimilitud, se usan frecuentemente. Por ejemplo,
el χ 2 de Pearson para independencia en tablas de contingencia es una
prueba escore, mientras el conocido estadístico G 2 para la misma situación
es una prueba del cociente de verosimilitud.

Para construir intervalos de confianza para parámetros se puede usar el

método de Wald o el método de verosimilitud. El primero es simplemente
βˆ ± zα s.e.βˆ , mientras que el último consiste en invertir la prueba de
2

hipótesis antes presentada.

Para la selección de covariables, en un proceso de modelación, es necesario

utilizar criterios de inclusión (o exclusión) basados en la función de
verosimilitud. Para definir un criterio de inclusión de una variable se debe
plantear un balance entre su “impacto” en la descripción del valor
esperado y la información importante que ese hecho necesita para el
proceso de estimación. Un estadístico apropiado para estos casos, para
comparar la bondad del ajuste del modelo y su grado de complejidad, es el
criterio de información de Akaike (Akaike, 1973), definido:
~
AIC p = −2 Lˆ ( p ) + 2 p = S p + 2 p − 2 L( n ) ,
~
Donde Lˆ ( p ) y L( n ) son las log-verosimilitudes evaluadas en βˆ ( p ) y
~
β ( n ) respectivamente y Sp es la deviance “escalada” (dividido por el
parámetro de escala estimado) para el modelo con p parámetros. Si el
modelo posee el parámetro de escala φ desconocido, previamente hay que
estimarlo. Este criterio extiende el criterio de cociente de máxima

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
verosimilitud para la situación de ajuste de muchos modelos con distinto
número de parámetros, decidiendo hasta cuándo se puede excluir o no
términos. Un valor bajo de AICp es considerado representativo de un
mejor ajuste y los modelos son seleccionados procurando obtener un
mínimo de AICp. Considerando dos modelos encajadas Mq y Mp, con p>q, se
tiene que AIC p AIC q = S p − S q − 2( p − q ) , así, suponiendo verdadero a Mq, se
prueba que E ( AIC p − AIC q ) = p − q + O(n −1 ).

2.3 Regresión de Poisson

Para modelar recuentos, la función de enlace identidad no es muy práctica

porque puede producir valores predichos negativos. La función logarítmica,
por otra parte, nos asegura que estos valores predichos van a ser
positivos. Por ejemplo,
log ( µ ( xi ) ) = α + β xi ; µ ( xi ) = exp(α + β xi ).

Podemos interpretar el efecto de la pendiente β en términos

multiplicativos: si x aumenta en 1 unidad, el promedio se multiplica por e β :

µ ( x + 1) = eα ( e β ) e β .
x

Similarmente, para variables regresoras cualitativos, como es el caso de

tener un modelo de análisis de la varianza con un fctor de tratamientos, la
diferencia entre los coeficientes asociados a dos tratamientos nos indica
por cuánto se multiplica la media al pasar de un tratamiento al otro. Para
el modelo log ( µi ) = α + τ i , al pasar, por ejemplo, del tratamiento 1 al
tratamiento 3, el promedio cambia: µ3 = µ1 exp (τ 3 − τ 1 ) . Esto quiere decir
que si τ 3 > τ 1 , entonces el factor es mayor que 1, y por lo tanto µ3 > µ1 . El
ajuste, la verificación del modelo y las pruebas de hipótesis se realizan de
la manera descripta para regresión logística.

3. Fenómeno de Superdispersión (variables discretas)

Dentro de la familia exponencial uniparámetrica, tanto la distribución

binomial como la Poisson, poseen la propiedad siguiente: la varianza queda
automáticamente definida a partir de la media. En efecto, en el caso
Poisson se tiene que var(Yi ) = µi . Sin embargo existen situaciones reales en
las que la varianza es más grande de lo que esperaríamos bajo el modelo.
Es más, McCullagh & Nelder (1989) afirman que en la práctica esto es lo

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
común y la expresión nominal, la excepción. Este fenómeno se llama
sobredispersión (superdispersión) y entre las causas más comunes
tenemos la heterogeneidad entre unidades de medición, o en el tiempo
como factor que indexa las respuestas. Cuando esto se presenta, hay dos
enfoques para la modelación de los datos:

• Proponer otra distribución más además de la elegida para la variable

aleatoria de respuesta, como si fuera una distribución compuesta,
de manera tal que el valor esperado sea considerado también como
variable (ej. binomial negativa o beta-binomial, en el caso de
variables binomiales, gama-Poisson, en el caso de conteos, etc.),
• Ajustar un modelo para la varianza mediante un factor, φ , y
estimarlo.

Para estimar φ , es necesario utilizar que la deviance o el χ 2 de Pearson se

aproximan (en muchos casos) con una distribución χ 2 . Por lo tanto, si el
modelo ajustara bien y no existiera sobredispersión, esperaríamos que
tanto la deviance como el χ 2 de Pearson estén cerca de los grados de
libertad (que son el valor esperado de una distribución χ 2 ). Ahora,
supongamos que el modelo sí ajusta bien (es decir, no hay falta de ajuste)
pero hay sobredispersión. Podemos demostrar que, asintóticamente, el
valor esperado de la deviance o del χ 2 de Pearson van a estar
multiplicados por φ . Esto hace que dos estimadores posibles de φ sean el
cociente entre la deviance y sus grados de libertad, o el cociente entre el
estadístico χ 2 de Pearson y sus grados de libertad. Hay que destacar que
esto será válido sólo si no hay falta de ajuste.

Este enfoque permite trabajar por analogía con σ 2 , el parámetro de

dispersión en la distribución normal. Para estimarlo usamos el cuadrado
medio residual (que es la deviance dividida por sus grados de libertad), lo
cual funcionará bien sólo en el caso de que no haya falta de ajuste. Los
estimadores del modelo se obtienen de la manera usual (como en la
distribución normal, los β̂ no dependen de σ 2 ) y la matriz de covarianza
de los mismos se multiplica por φˆ (recordemos que la covarianza de los β̂
es σ 2 ( X ' X ) ). Luego, con todo, las pruebas de hipótesis (Wald, Máxima
−1

Verosimilitud y Escore) quedan divididas por φˆ (es decir, sean más

conservadoras) y los errores estándar asintóticos quedan multiplicados

1
Basada en la bibliografía sugerida en la planificación del curso. 2009
por φˆ . Por analogía al modelo clásico normal y el análisis de la varianza de
Fisher, se ha propuesto usar la distribución F en vez de la χ 2 , ya que φ es
estimado.

4. Regresión Poisson con offset

Existen muchas aplicaciones en las que no nos interesa modelar los
recuentos propiamente dichos sino una tasa, densidad, promedio, etc. Por
ejemplo, supongamos que contamos, en 50 plantas tratadas y en 50 plantas
consideradas como controles, la cantidad de insectos (Y) y la cantidad de
hojas (h), y nos interesa modelar el promedio de insectos por hoja. Un
µi
modelo que podemos usar es log = α + τ i ,o log µi = log hij + α + τ i . El primer
hij
término se llama offset. Las características del MLG no cambian, es decir
la estructura de predictor lineal, función de enlace y distribución de la
respuesta. Observemos que simplemente es como agregarle una nueva
variable independiente al modelo, cuyo coeficiente de regresión es 1).

5. Regresión Binomial Negativa

Si bien la distribución binomial negativa no es una familia exponencial para
k desconocido, la cantidad de comportamientos para datos biológicos que
requieren de esta función de distribución de probabilidad hace que se
hayan desarrollado métodos que permiten usar las técnicas de modelos
lineales generalizados para datos binomiales negativos. Esta distribución
es descripta por la siguiente expresión:

P( y) =
(
Γ y+ 1
k ) (kµ )
k

y = 0,1, 2,...; µ > 0; k > 0.

(
Γ ( y + 1) Γ 1
k ) (1 + k µ )
y+ 1
k

El parámetro k se refiere comúnmente como el índice de “agregación”. La

µ2
función de varianza es Var(Y ) = µ + , con lo que vemos que ésta es una
k
distribución alternativa para modelar sobredispersión en recuentos
Poisson. Si k es conocido, se puede modelar directamente ya que es una
familia exponencial. Si k es desconocido, la estrategia es estimarlo
conjuntamente con los parámetros del modelo mediante máxima
verosimilitud (análogo a estimar σ 2 junto con los otros parámetros en
regresión normal). La función de enlace usada típicamente es la
logarítmica.

1
Basada en la bibliografía sugerida en la planificación del curso. 2009