ESTADÍSTICA ESPACIAL EN EPIDEMIOLOGÍA Y MEDIO AMBIENTE

DOCTORADO EN ESTADÍSTICA E I.O.

Universitat de València (Estudi General)
Burjassot, Primavera de 2004

MODELOS LINEALES GENERALIZADOS

Antonio López
Dep. d’Estadı́stica i Investigació Operativa
Universitat de València (Estudi General)
[Link]@[Link]

1
GUIÓN:
1 Introducción
Motivación. Mortalidad por cáncer de próstata en Valencia. Regresión
lineal.

2 Modelo Lineal Generalizado (GLM)

Definición. Componentes. Función vı́nculo. Modelos de datos
continuos. Modelos de datos discretos. Parámetro de dispersión.
Sobredispersión.

3 Estimación de un GLM
Máxima Verosimilitud. Método Scoring de Fisher. Estimación del
parámetro de dispersión.

4 Selección del mejor modelo

Desviación. Modelos encajados.

2
5 Análisis de residuos
Residuos de Pearson. Residuos de desviación. Residuos por exclusión.

6 Extensiones de los GLM

Quasi-verosimilitud. Otras extensiones.

7 Ejemplo
Mortalidad por cáncer de próstata en Valencia.

Bibliografı́a

3
MORTALIDAD POR CÁNCER DE PRÓSTATA.
VALENCIA 1975-1980

0 6 a 10
1a5 >10

número de defunciones acumuladas para el perı́odo

4
TASAS DE MORTALIDAD POR CÁNCER DE
PRÓSTATA. VALENCIA 1975-80

0 [5,10[
]0,5[ [10,150[

tasas por 10000 habitantes

5
CONCENTRACIÓN DE NITRATOS EN AGUAS
POTABLES

[ 0,10[ [30,90[
[10,30[ [90,300[

Concentración de nitratos en mg/litro, Llopis (1985)

6
 MODELOS DE REGRESIÓN

yi = f (xi ) + εi i = 1, . . . , n indep.

atributo en componente componente

= +
estudio sistemática errática

Esperanza[yi ] = f (x)
Varianza[yi ] = Varianza[εi ]

Regresión lineal simple:

f (xi ) = β0 + β1 xi
Ej.: mortalidadi =tasa×poblacióni +εi
i:= ı́ndice de municipio

7
 MODELOS DE REGRESIÓN

yi = f (xi ) + εi i = 1, . . . , n indep.

atributo en componente componente

= +
estudio sistemática errática

Esperanza[yi ] = f (x)
Varianza[yi ] = Varianza[εi ]

Regresión lineal múltiple:

f (xi ) = β0 + β1 x1i + . . . + βk xki
P
Ej.: mortalidadi = j (tasaj ×poblaciónji ) +εi
i:= ı́ndice de municipio
j:= ı́ndice de grupo de edad

8
C. DE PRÓSTATA EN VALENCIA (cont.)

NOMBRE DEL cáncer conc. edad:

MUNICIPIO habit. prostata nitr. % ≥ 40

Ademuz 1545 1 11 59.0

Ador 1256 5 16 49.4
Adzaneta de Albaida 1364 0 18 42.7
Agullent 2016 0 8 35.8
Alaquàs 23728 5 78 32.4
Albaida 5573 3 8 38.7
Albal 8139 4 17 36.0
Albalat de la Ribera 3594 2 76 42.2
Albalat dels Sorells 567 8 60 41.4
Albalat dels Tarongers 3657 0 32 53.4
Alberique 8971 1 28 42.9
Alborache 821 0 12 43.9
Alboraya 10786 4 42 39.2
Albuixech 3005 0 66 47.0
... ... ... ... ...

Datos de nitratos extraı́dos de (Llopis, 1985)

9
 REGRESIÓN LINEAL SIMPLE
σ2 )
yi ∼ N(β0 + β1 xi , |{z} i = 1, . . . , n
| {z }
media varianza independ.

A partir de los n datos se obtienen las estimaciones:

P
i (y
Pi − ȳ)(xi − x̄)
β̂1 =
n i (xi − x̄)2
β̂0 = ȳ − β̂1 x̄

y se contrasta la hipótesis H0 : β1 = 0

10
 C. DE PRÓSTATA EN VALENCIA (cont.)

casosi = β0 + β1 × nitratosi + εi

Ajuste
σ̂ 2 : 16.88 con 261 grados de libertad
Multiple R2 : 0.00007167
Estadı́stico F : 0.01871 con 1 y 261 gr. libertad,
(p-valor=0.8913) resultado no significativo

Coeficientes Valor [Link]. estad.t p-valor

Intercept. 3.0091 1.4558 2.0669 0.0397
nitratos 0.0032 0.0234 0.1368 0.8913

11
 C. DE PRÓSTATA EN VALENCIA (cont.)

Residuos:
Min 1Q Median 3Q Max
-3.534 -3.06 -2.102 -0.3581 265.9

12
 C. DE PRÓSTATA EN VALENCIA (cont.)

ajuste de mortalidad vs. nitratos

o
200
mortalidad

100
50

o oo ooo o
o o o
ooo
ooooooooooo
o o o
oooooooo
oo
oo o o o
oooooooooooooooooo
oooooooooooooo ooo
ooo oo ooooo
ooooooooooooo o o
oooo o o oooo oo oooo o ooo o o
0

0 50 100 150 200 250

nitratos

13
 C. DE PRÓSTATA EN VALENCIA (cont.)

residuos del ajuste eliminando municipio de Valencia

o
20

o
o
15

o
o
residuos

o
10

o
o o o o
o o
o o
o o o
5

o
oo o o o ooo ooo o o
o o o o o o
o ooooo o o oo o oo o o o oo oo o oo o ooo
o oo o
0

oo o oooo oooo o oooo o o oo o o

o oo o oo o ooo o o oooooooo o oo oooooooooooooo o o oo oo ooooooooo o oooo
o oooooooooo oooo ooooo oo ooo o oooooooo ooo oooooo oooooo oooooooo oooooooooooooooo ooooooo oooooooo
o o o oo o
o o
oo
o o
o oo o o o o ooo oo o o

0 50 100 150 200 250

municipios

14
 REGRESION LINEAL MULTIPLE
Aunque estemos interesados en un sólo factor de riesgo, conviene
incluir en el modelo todos aquellos cuya influencia se sospecha. Con
ello evitaremos estimaciones sesgadas del factor de interés y
conclusiones posiblemente equivocadas.
Ajustando
yi = α0 + α1 x1i + εi
y el modelo con factores de riesgo x2 , . . . , xk adicionales
yi = β0 + β1 x1i + . . . + βk xki + εi
en general se obtiene α̂1 6= β̂1
La tabla de ANOVA permite valorar la hipótesis de que todos los
coeficientes de regresión son nulos a la vez, ası́ como la colección de
hipótesis de que cada uno de ellos es irrelevante una vez considerados
los restantes.
15
 C. DE PRÓSTATA EN VALENCIA (cont.)

casosi = β0 + β1 × poblacióni
+ β2 × envejecimientoi
+ β3 × nitratosi + εi

Ajuste
σ̂ 2 : 2.058 con 259 grados de libertad
Multiple R2 : 0.9852
Estadı́stico F : 5765 con 3 y 259 grad. de lib,
(p-valor < 0,00005) muy significativo

16
 C. DE PRÓSTATA EN VALENCIA (cont.)
Coefs. Estim. StdErr. estad.t p-valor
Intercept. -1.5897 0.8538 -1.8618 0.0638
población 0.0004 0.0000 130.60 0.0000
envejec. 3.5904 1.7374 2.0665 0.0398
nitratos 0.0061 0.0030 2.0291 0.0435

Residuos
Min 1Q Median 3Q Max
-8.765 -0.7577 -0.3334 0.5301 11.65
Incremento R2 ≈ 0.9851

17
 DIAGNÓSTICO DEL MODELO
ajuste de mortalidad vs. pobl., envej. y nitratos

o
10

oo o
o
o o o o o oo
5

o
residuos

o o o oo o o
o o o o o o
oo o o oooo oo oooooooo o oooo o ooo ooo oooo ooooo ooo oooo oo o
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo oooooooooooooooooooooooooooooooooo oooo oooooooo ooooooo ooooooooooooooooooooooooooooo
0

o oo o o oooooo oooo
ooo o o oo o o oo o o o ooo o o
o o o o o
o oo o o
o
-5

o o
o

0 50 100 150 200 250

municipio

18
 DIAGNÓSTICO DEL MODELO
mortalidad ajustada para demografia vs. nitratos

o
10

o oo
mortalidad ajustada

o o o
oo o o
5

o o
o o
oo oo ooooo oo o o o ooo o o o o
ooooooooooooo ooo oo o oo o ooo o oo o
oooooo
o
o
oooo
oooo oooo ooooooooo o oooooo ooo
oooooo o ooo oo o
0

oo
o
oo
oo o
oo
o
oo
ooo
o
oo
oo
oo
oo
oo
ooo oooooo o oo
o ooo oo ooooooooooo oo o oo o o
ooo o oo o oo o o o o
o o o o o oo o
o o
-5

o o
o

0 50 100 150 200 250

nitratos

19
 Introducción
Modelo Lineal General: datos independientes, y1 , y2 , . . . , yn ,
normalmente distribuidos.

yi ∼ N(β0 + β1 x1i + . . . + βp xpi , σ 2 )

predictor lineal β 0 xi
varianza constante

E[y] = Xβ, V[y] = σ 2 I

Modelo Lineal Generalizado: datos independientes de una

distribución de la familia exponencial (binomial, Poisson,
gamma, . . .).
modeliza E[y] como una función no lineal de Xβ.

20
 Introducción

Análisis de un GLM:
cálculo del estimador máximo verosı́mil
comparación de modelos encajados
valoración del ajuste del modelo a los datos

21
 Definición de GLM
Conjunto de variables aleatorias independientes y1 , y2 , . . . , yn con
función de densidad, o función de probabilidad, que puede escribirse
como:
yi θi − b(θi )
p(yi | θi , φ) = exp{ + c(yi , φ)}
ai (φ)
donde:
θi es el parámetro natural o canónico
φ es un parámetro adicional de escala o dispersión
ai (·), b(·) y c(·) son funciones especı́ficas
Si φ es conocido este es un modelo de la familia exponencial lineal
Si φ es desconocido es un modelo de dispersión exponencial

22
 Definición de GLM

Bibliografı́a general:
Nelder y Wedderburn (1972)
McCullagh y Nelder (1989)
Fahrmeir y Tutz (1994)
Garthwaite et al. (1995)

23
 Componentes del GLM
Queremos modelizar µi = E[yi ] en términos del predictor lineal β 0 xi
formado con un conjunto de p covariables

β 0 xi = β0 + β1 x1i + . . . + βp xpi

Componentes:
1 Conjunto de n variables respuesta independientes, de una
distribución de la familia exponencial
2 Un vector de parámetros β y una matriz del modelo X,
determinando el predictor lineal de cada variable β 0 xi
3 Una función vı́nculo monótona y diferenciable que define la
relación entre µi y su predictor lineal

g(µi ) = β 0 xi

24
 Función vı́nculo
Permite modelizar distintas relaciones entre µ y el predictor lineal.
Vı́nculo natural o canónico:
Aquel que es igual a la función que define el parámetro natural o
canónico de esa distribución. Por tanto, θ = β 0 x

25
 Función vı́nculo
Vı́nculos más usuales:
π
¦ logit log 1−π
−1
¦ probit Φ (π)
¦ complementario
log-log log[− log(1 − π)]
¦ identidad µ
¦ inverso −1/µ
¦ logaritmo log µ
√
¦ raiz cuadrada µ
Elección del vı́nculo: depende de la familia de distribuciones, del
tipo de respuestas y de la aplicación.

26
 Modelos de datos continuos
Normal:

Distribución N(µ, σ 2 )

E[y] = µ

vı́nculo g(µ) = µ (identidad)

b(θ) = θ2 /2

a(φ) = σ2

Otros vı́nculos: logaritmo

raiz cuadrada

27
 Modelos de datos continuos
Gamma:

Distribución Gamma(λ, ν)

λ
E[y] = ν

vı́nculo g(µ) = − µ1 = − λν (inverso)

b(θ) = − log(−θ)

1
a(φ) = λ

Otros vı́nculos: identidad

logaritmo

28
 Modelos de datos discretos
Binomial:

Distribución Bi(n, π)

E[y] = nπ

µ π
vı́nculo g(µ) = log n−µ = log 1−π (logit)

b(θ) = n log(1 + eθ )

a(φ) = 1

Otros vı́nculos: probit

complementario log-log

29
 Modelos de datos discretos
Poisson:

Distribución Po(λ)

E[y] = λ

vı́nculo g(λ) = log λ (logaritmo)

b(θ) = eθ

a(φ) = 1

Otros vı́nculos: identidad

raiz cuadrada

30
 Parámetro de dispersión
Con frecuencia, el término ai (φ) es de la forma φ/ωi , donde ωi es un
peso.
Si los datos no son agrupados, ωi = 1

Si las variables respuesta expresan promedios, ωi = ni

Si son la suma de ni respuestas individuales, ωi = 1/ni

31
 Sobredispersión
Fenómeno que ocurre en aplicaciones con distribuciones con varianza
poco flexible, como Binomial y Poisson.
Al añadir un parámetro de dispersión φ, se modifica la varianza

V[y] = a(φ)b00 (θ)

Puede representar una heterogeneidad no observada o una correlación

positiva entre respuestas individuales.
También se denomina extravarianza.

32
 Máxima verosimilitud
El logaritmo de la verosimilitud de θ para las observaciones y es
n
X n
X
yi θi − b(θi )
l(θ | y) = + c(yi , φ)
i=1
ai (φ) i=1

Nuestro principal interés es la estimación de β. El estimador

máximo verosı́mil de cada βj anula la derivada de l
n
X (yi − µi )xij
∂l
= 0 (µ )
∂βj i=1
V[yi ]g i

33
 Máxima verosimilitud
En general, estas ecuaciones de estimación no se pueden resolver
directamente. Su solución puede aproximarse por procedimientos
iterativos, empleando la esperanza de las segundas derivadas
· ¸ X n
∂2l xij xik
E = 0 (µ )2
∂βj ∂βk i=1
V[y i ]g i

34
 Método Scoring de Fisher
Algoritmo de Newton-Raphson:
Procedimiento iterativo a partir de una estimación inicial β 0 :
β r+1 = β r − [Dβ2 l(β r )]−1 Dβ l(β r )
donde Dβ l(β r ) es el vector de primeras derivadas de l, y Dβ2 l(β r )
la matriz de segundas derivadas, evaluadas en β r .
Método Scoring de Fisher:
Consiste en sustituir Dβ2 l(β r ) por su valor esperado.
· 2
¸ X n
∂ l xij xik
E = 0 (µ )2
∂βj ∂βk i=1
V[y i ]g i

Equivale a resolver iterativamente un problema de mı́nimos

cuadrados ponderados (Jorgensen, 1983).
La sucesión {β r } converge al estimador máximo verosı́mil de β.

35
 Estimación del
parámetro de dispersión
Si φ no es conocido, es necesario usar una estimación para el cálculo
de V[yi ] en el procedimiento anterior.
Cuando ai (φ) = φ/ωi , la expresión de la varianza

V[yi ] = ai (φ)b00 (θi )

proporciona un estimador consistente de φ a partir de una

estimación de β
n
X ωi (yi − µ̂i )2
1
φ̂ =
n − p − 1 i=1 b00 (θ̂i )

36
 Estimación del
parámetro de dispersión
Para la normal, el estimador de la varianza del modelo de regresión
lineal es la suma de cuadrados residual
X n
2 1
σ̂ = (yi − µ̂i )2
n − p − 1 i=1

37
 Desviación
Determinaremos la adecuación del modelo comparándolo con el
modelo saturado.
El modelo saturado tiene la misma forma que el ajustado, pero
con tantos parámetros como observaciones.
Desviación escalada: obtenida con el estadı́stico cociente de
verosimilitudes

S = −2[l(β̂ | y, φ) − l(β̃ | y, φ)]

con β̃ el EMV del modelo saturado.

38
 Desviación
En términos del parámetro natural es
n
X yi (θ˜i − θˆi ) − b(θ˜i ) + b(θˆi )
S=2
i=1
ai (φ)

Cuando φ es conocido, la desviación escalada mide cuánto se desvı́a

el modelo de los datos.
Distribución aproximada:
Si el modelo se ajusta bien a los datos

S ∼ χ2 (n − p − 1)

39
 Desviación
Desviación (no escalada):
Se define por
D(y, µ̂) = φS

Si ai (φ) = φ/ωi , equivale a

n
X
2 ωi [yi (θ˜i − θˆi ) − b(θ˜i ) + b(θˆi )]
i=1

Descomposición de la desviación:
La desviación es la suma de las discrepancias para cada uno de
los datos
n
X
D(y, µ̂) = di (yi , µ̂i )
i=1

40
 Desviación
Estimación de φ:
La desviación de un modelo razonable con q parámetros permite
estimar φ mediante
φ̂ = D/(n − q)
debido a que la esperanza aproximada de S es igual a n − q, los
grados de libertad de la distribución χ2

41
 Modelos encajados
La desviación es útil para comparar el ajuste de dos modelos
encajados.
Un modelo M1 con q1 parámetros está encajado en otro M2 con q2
parámetros (q1 < q2 ) si son de la misma forma y las covariables
de M1 están contenidas en las de M2 .
La necesidad de los q2 − q1 parámetros adicionales se contrasta con
un test χ2 . Si D1 y D2 son las desviaciones de dos modelos
encajados con buen ajuste,

(D1 − D2 )/φ ∼ χ2 (q2 − q1 )

42
 Modelos encajados
Si φ tiene que ser estimado, puede hacerse el contraste con un test
F, usando
(D1 − D2 )(n − q2 )
∼ F(q2 − q1 , n − q2 )
(q2 − q1 )D2

43
 Análisis de residuos
El residuo de cada dato mide la discrepancia entre el valor observado
y el pronosticado por el modelo.
Residuos de Pearson: Generalización inmediata de los residuos
habituales para datos normales
yi − µ̂i
riP =q
b00 (θ̂i )

Residuos de desviación: Es la contribución de esa observación a

la desviación escalada
D
p
ri = signo(yi − µ̂i ) di /φ

Residuos por exclusión: Es el residuo de ese punto para el modelo

ajustado al excluir esa observación. Pueden calcularse residuos
por exclusión de Pearson y de desviación.

44
 Quasi-verosimilitud
A veces no se conoce la forma de la distribución de las variables
respuesta, pero se dispone de la esperanza en función de β
E[yi ] = µi (β)
y la fórmula de la varianza en su relación con la esperanza
V[yi ] = φV(µi )
Estimador por quasi-verosimilitud
Es la solución de
D0 W (y − µ(β)) = 0
∂µi
donde el elemento (i, j) de D es ∂βj y W es la matriz diagonal
con elementos V(µi )−1 .
Quasi-desviación
Como la desviación, sustituyendo por la quasi-verosimilitud.

45
 Otras extensiones
Modelos de regresión no lineal
Empleando un predictor no lineal en los parámetros β.
Modelos de regresión general
Utilizando distribuciones que no son de la familia exponencial.
Modelos de regresión multivariante
La variable respuesta es un vector, introduciendo los GLM
multivariantes
(Fahrmeir y Tutz, 1994).
O las respuestas no son independientes, como en el caso espacial,
llevando a los modelos autoregresivos y a los jerárquicos.

46
 Mortalidad por cáncer de próstata en Valencia

Estimación del modelo.

Parámetros estimados
MODELO β0 β1 β2
tasas const. -7.172
edad -9.925 5.208
nitratos -7.876 1.23e-3
edad y nit. -10.152 5.539 2.09e-3

47
 Mortalidad por cáncer de próstata en Valencia

Diferencias entre las desviaciones de los modelos encajados.

const. edad nit. comp.
tasas const. 849.8
edad 488* 361.8
nitratos 443* — 406.8
edad y nit. 495.9* 7.9* 52.9* 353.9

Todas significativas con α = 0,01.

48
 Bibliografı́a
Fahrmeir, L. y Tutz, G. (1994). Multivariate statistical modelling based
on generalized linear models. Springer-Verlag, New York.
Ferrándiz, J., López, A., Llopis, A., Morales, M., y Tejerizo, M. L.
(1995). Spatial interaction between neighbouring counties: cancer
mortality data in Valencia, (Spain). Biometrics, 51(2):665–678.
Garthwaite, P. H., Jolliffe, I. T. y Jones, B. (1995). Statistical Inference.
Prentice Hall, London.
Jorgensen, B. (1983). Maximum likelihood estimation and large-sample
inference for generalized linear and nonlinear regression models.
Biometrika, 70:19–28.
McCullagh, P. y Nelder, J.A. (1989). Generalized linear models, second
edition. Chapman and Hall, London.
Nelder, J.A. y Wedderburn, R.W.M. (1972). Generalized linear models.
Journal of the Royal Statistical Society, series A, 135:370–384.

49