Capı́tulo 5

Otros modelos glm: multinomial,

ordinal y poisson

1. Regresión Multinomial
Hasta ahora nos hemos centrado en el caso en el que la variable respuesta era dicotómica.
Ahora nos centramos en el caso en el que la variable de interés tiene más de dos categorı́as,
por simplicidad, vamos a ilustrar la metodologı́a para el caso de tres categorı́as, ya que la
generalización a más de tres es inmediata.

Supongamos que codificamos las tres categorı́as de la variable respuesta como 0, 1 y 2.

En el caso de regresión logı́stica, el logit es
 
P r(Y = 1)
log
P r(Y = 0)

Ahora el modelo necesita dos funciones logit ya que tenemos tres categorı́as, y necesitamos
decidir que categorı́as queremos comparar. Lo más general es utilizar Y = 0 como referencia
y formar logits comparándola con Y = 1 y Y = 2. Supongamos que tenemos k variables
explicativas, entonces:
 
P r(Y = 1)
log = β10 + β11 X1 + β12 X2 + . . . + β1k Xk
P r(Y = 0)
 
P r(Y = 2)
log = β20 + β21 X1 + β22 X2 + . . . + β2k Xk
P r(Y = 0)

1
con lo cual ahora tenemos el doble de coeficientes que en el caso de regresión logı́stica. Las
probabilidades se calcuları́an como:
1
P r(Y = 0|X) =
1 + eg1 (X)+g2 (X)
eg1 (X)
P r(Y = 1|X) =
e1+g1 (X)+g2 (X)
eg2 (X)
P r(Y = 2|X) =
e1+g1 (X)+g2 (X)
g1 (X) = β10 + β11 X1 + β12 X2 + . . . + β1k Xk
g2 (X) = β20 + β21 X1 + β22 X2 + . . . + β2k Xk

La estimación de los parámetros se hace nuevamente mediante el método de máxima verosimi-

litud.

Los datos que vamos a utilzar en este caso provienen de un estudio cuyo objetivo es es-
tudiar los factores que influyen en el conocimiento, actitud y comportamiento de mujeres
hacia las mamografı́as. Son datos sobre 412 mujeres y las variables son:

Variable Descripción Valores

me Frecuencia del chequeo 0=Nunca
1=El último año
2= Más de un año
symp La mamografı́a no es 1=Muy de acuerdo
necesaria a menos que 2=De acuerdo
tengas algún sı́ntoma 3=En desacuerdo
4=Muy en desacuerdo
pb Beneficio percibido 5-20
hist Madre y hermana con 0=No
historia de cancer de mama 1=Sı́
bse ¿Te han enseñado a explorarte? 0=No
1=Sı́
detc ¿Es probable que una 0=No
mamografı́a detecte un 1= Algo
nuevo caso de cancer? 2= Muy

mamexp=read.table( "mamexp.txt",header=TRUE)
names(mamexp)
mamexp$symp=factor(mamexp$symp)
mamexp$hist=factor(mamexp$hist)
mamexp$bse=factor(mamexp$bse)
mamexp$dect=factor(mamexp$dect)

attach(mamexp)
table(me)

2
me
0 1 2
234 104 74

1.1. El procedimiento multinom() en R

La sintaxis del comando multinom es similar a la de los comandos para regresión logı́stica.
Para utilizar esta función necesitamos cargar una nueva librerı́a

library(nnet)
multi1=multinom(me~hist)
multi1
Coefficients:
(Intercept) hist1
1 -0.9509794 1.256531
2 -1.2504803 1.009384

logit(Último año/Nunca|hist) = −0,95 + 1,25 × hist

logit(Más de un año/Nunca|hist) = −1,25 + 1,01 × hist

Por defecto, el nivel estamos tomando como referencia es el primer nivel. El output es similar
al de regresión logı́stica, sin embargo, tenemos dos grupos de parámetros, uno para cada logit.
En este caso, no tenemos un p-valor para los coeficientes, si queremos saber si la variable
explicativa es significativa, utilizamos el test de la razón de verosimilitud

multi0=multinom(me~1)
multi0$deviance-multi1$deviance
multi1$edf-multi0$edf
1-pchisq(12.858,2)

1.2. Interpretación y significación de los parámetros

El OR en este caso serı́a
P r(Y = j|X = a)/P r(Y = 0|X = a)
ORj (a, b) = j = 1, 2...
P r(Y = j|X = b)/P r(Y = 0|X = b)

Donde j indica el nivel de la variable respuesta con el que estamos comparando y a y b

son los valores que toman las variable explicativas. Si sólo tenemos una variable explicativa
dicotómica que toma valores 0 y 1, entonces tendrı́amos Rj (1, 0). En el ejemplo:

table(me,hist)
hist

3
me 0 1
0 220 14
1 85 19
2 63 11

d 1 = 19 × 220 = 3,51 = exp(1,25)

OR
85 × 14
11 × 220
OR
d2 = = 2,74 = exp(1,009)
63 × 14
La interpretación de los parámetros, los OR es similar al caso de regresión logı́stica. La
interpretación del efecto de la historia familiar en la frecuencia de mamografı́as es la siguiente:

Mujeres con historia familiar de cancer de pecho tienen un odds de haberse hecho una
mamografı́a en el último año que es 3.51 veces mayor que el de las mujeres sin historia
familiar.

Mujeres con historia familiar de cancer de pecho tienen un odds de haberse hecho
una mamografı́a hace más de un año que es 2.7 veces mayor que el de las mujeres sin
historia familiar.

En otras palabras, tener un familiar directo con cancer de pecho es un factor significativo en
el uso de mamografı́as.

Si en vez de hist introducimos dect (la opinión de las mujeres sobre la habilidad de la
mamografı́a para detectar el cancer):

multi2=multinom(me~dect)
summary(multi2)
exp(coef(multi2))

Si nos fijamos en los errores estándar, vemos que todos los intervalos de confianza con-
tendrı́an al 0 excepto uno, ese OR se interpretarı́a de la siguiente forma: El odds de hacerse

4
una mamografı́a dentro del último año entre las mujeres que piensan que es muy probable
que la mamografı́a detecte un nuevo caso de cáncer es 8.22 veces mayor que el odds entre
las mujeres que piensan que es no es probable que la mamografı́a lo detecte.

Si la variable explicativa es continua, aparece un parámetro para cada logit y representa

el OR para un cambio en una unidad de la variable continua.

1.3. Selección de variables

Los pasos a seguir serı́an los mismos que para regresión logı́stica, comenzarı́amos viendo
si cada una de las variables por separado son significativas. Si tenemos los modelos
multi3=multinom(me~symp)
multi4=multinom(me~bse)
multi5=multinom(me~pb)
¿Cuáles son significativas?
multi6=multinom( me~hist+dect+symp+bse+pb)
summary(multi6)
Call:
multinom(formula = me ~ hist + dect + symp + bse + pb)

Coefficients:
(Intercept) hist1 dect1 dect2 symp2 symp3 symp4
1 -2.9990607 1.366274 0.01700706 0.9041717 0.1101408 1.9249158 2.457230
2 -0.9860258 1.065446 -0.92448094 -0.6906153 -0.2901346 0.8172916 1.132224
bse1 pb
1 1.291772 -0.2194421
2 1.052183 -0.1482079

Std. Errors:
(Intercept) hist1 dect1 dect2 symp2 symp3 symp4
1 1.539293 0.4375239 1.1619394 1.1268643 0.9228324 0.7776651 0.775400
2 1.111829 0.4593986 0.7137332 0.6871025 0.6440622 0.5397872 0.547666
bse1 pb
1 0.5299080 0.07551477
2 0.5149934 0.07636886

Residual Deviance: 693.9019

AIC: 729.902
Los dos coeficientes estimados para symp que estiman el log-odds de las que están de acuerdo
frente al grupo de referencia que son las que están muy de acuerdo no parecen ser signifi-
cativos, por lo que podrı́amos simplificar la codificación de la variable agrupando los dos
primeros niveles.
Vemos además que ninguno de los cuatro coeficientes de dect parece ser significativo, por

5
lo que decidimos ajustar el modelo sin esta variable y nos fijamos en dos cosas: 1) si los
parámetros de las demás variables cambian (con lo cual estamos comprobando si es una
variable de confusión) y 2) qué nos dice el test de la razón de verosimilitud.
multi6=multinom( me~hist+dect+symp+bse+pb)
multi7=multinom( me~hist+symp+bse+pb)
coef(multi6)
coef(multi7)
multi7$deviance-multi6$deviance
multi6$edf-multi7$edf
1-pchisq(8.44,4)

En cuanto a la variable pb vamos a codificarla para poder comprobar si la relación es

lineal:
benef=pb
benef[benef<=5]=0
benef[benef>5 & benef<=7]=1
benef[benef>7 & benef<=9]=2
benef[benef>9]=3
benef=factor(benef)
multi8=multinom( me~hist+symp+bse+benef)
hacemos un gráfico de los coeficientes para los dos logits:
coefi=coef(multi8)
coefi=coefi[,7:9]
matplot(t(coefi),type="l",xlab="Niveles",ylab="Coeficientes")
El gráfico confirma que la relación es lineal, con lo cual tendrı́amos la opción de dejar la
variable continua o utilizar la categorizada.

6
 −0.4
−0.6
Coeficientes

−0.8
−1.0
−1.2
−1.4

1.0 1.5 2.0 2.5 3.0

Niveles

Ahora tenemos que comprobar si las interacciones son significativas, introduciendolas una a
una en el modelo. En este caso, ninguna de las diez posibles intercacciones es significtiva.

Los OR estimados:
exp(coef(multi8)
(Intercept) hist1 symp2 symp3 symp4 bse1 benef1
1 0.03145177 3.862051 1.1079680 7.822312 13.657506 3.542514 0.7474058
2 0.09006473 3.005730 0.7513073 2.355203 3.291352 2.710528 0.7086432
benef2 benef3
1 0.4828098 0.2330417
2 0.5961067 0.4704700
Los OR estimados muestran que hay un aumento en el odds de hacerse una mamografı́a
para ambas categorı́as de frecuencia frente a la categorı́a “nunca” (excepto para la variable
beneficio, ya que valores más altos de esta variable corresponden a una menor creencia en
el beneficio de este tipo de diagnóstico). Las estimaciones de los OR son mayores para -
eciente”frente a ”nunca”, lo que indica que los dos grupos perciben de forma diferente el
valor de la mamografı́a.
En el caso de la variable sı́ntoma, las mujeres que no están de acuerdo con la frase:
”La mamografı́a no es necesaria a menos que aparezca algún sı́ntoma”tienen un odds
13.65 veces mayor de haberse hecho una mamografı́a recientemente, y un odds 3.2 veces
mayor de habérsela hecho no tan recientemente, comparado con mujeres que no están
en desacuerdo con esa frase.
Las mujeres con historia familiar de cáncer tienen un odds 3.8 veces mayor de haber
hecho uso de la mamografı́a recientemente y 3 veces mayor en el caso no tan reciente
comparado con mujeres que no tienen un familiar próximo con cáncer.
Haber aprendido a explorarse el pecho es un factor significtivo a la hora de haberse
hecho una mamografı́a en el último año, pero no lo es a la hora de hacerlo no tan
recientemente (ver los errores estándar).

7
 La variable benef tiene un OR < 1 ya que valores más altos indica creencia de menor
beneficio, además, cuanto más aumenta la desconfianza más disminuye el odds de
someterse a una mamografı́a ya sea reciente o no recientemente.

Como ya mencionamos anteriormente, cuando utilizamos un modelo de regresión multino-

mial, la interpretación de los coeficientes es más complicada, ya que tenemos más de un OR
para cada variable. Sin embargo, utilizar una respuesta multinomial da una información más
completa sobre el proceso que estamos estudiando. Por ejemplo, si hubiéramos codificado
la variable respuesta como alguna vez frente a nunca habrı́amos perdido la información que
indica que los odds son mayores para el uso frecuente.
Desde el punto de vista estadı́stico, no es conveniente pasar a una variable dicotómica a
menos que los coeficientes en los distintos logits no se puedan considerar significativamente
distintos.

2. Regresión para datos ordinales

¿Cuándo unos datos son ordinales? Son datos que toman valores que tienen un orden
natural, para los cuales los intervalos entre valores no tiene por qué tener un significado. Un
ejemplo de datos ordinales: estado de salud: Excelente, muy bueno, bueno, regular, malo;
encuestas de satisfacción: Mucho, Poco, Nada, etc.

Vamos a trabajar con datos de un estudio cuyo interés es buscar los factores que influ-
yen en un graduado a la hora de decidir va a realizar estudios de postgrado o no. A 400
estudiantes de grado se les preguntó si era: nada probable, algo problable, o muy probable
que realizaran estudios de postgrado. Se recogieron datos sobre el nivel de educación de sus
padres (0=bajo, 1=alto), si el centro en el que estaban estudiando era público o privado, y
su nota de expediente académico.
Cuando nos encontramos con este tipo de datos podemos tomar varias decisiones:

1. Recodificar los datos, pasar a una variable dicotómica y utilizar los modelos que vimos
en el capı́tulo 2.

2. Analizar cada unas de las categorı́as por separado:Muy probable/Algo probable/Nada

probable

Estos modelos son fáciles de interpretar, pero son estadı́sticamente ineficientes ya que no
utilizan toda la información que proporcionan los datos. La regresión logı́stistica ordinal
combina los modelos anteriores en uno solo.

2.1. Modelo de odds proporcionales

Lo que este modelo hace es trasformar la escala ordinal a varios puntos de corte binarios,
el número de puntos de corte es siempre menor que el número de categrorı́as. Por ejemplo, si
hay tres categorı́as, habrı́a dos puntos de corte. Es lógico pensar estos puntos de corte como
umbrales que necesitamos cruzar para pasar una categorı́a a la siguiente (más alta) categorı́a.
Los dos umbrales en nuestro caso nos harı́a pasar de nada probable a algo probable, y de

8
algo probable a muy probable. En este modelo no estimamos la probabilidad de algo, sino
la probabilidad de observar un resultado o algo menor.
En general, supongamos que la variable respuesta Y toma valores de 1, . . . , K (en el
ejemplo, de 1 a 3), definimos

p̃k = P t(Y ≤ k|X), para k = 1 . . . K − 1

Las pk son las probabilidades acumuladas y k los puntos de corte.

En el modelo de odds-proporcionales se comparan la probabilidad de obtener una respuesta

<= k con la de tener una respuesta > k

 
p1
logit(p̃1 ) = log = α1 − β1 X1 + . . . − βm Xm
1 − p1
 
p1 + p2
logit(p̃2 ) = log = α2 − β1 X1 + . . . − βm Xm
1 − p1 − p2
 
p1 + p2 + . . . + p k
logit(p̃k ) = log = αk − β1 X1 + . . . − βm Xm
1 − p1 − p2 − . . . − pk
1 = p1 + p2 + . . . + pk + pk+1

El modelo asume que los coeficientes que acompañan a las covariables son iguales, y la única
diferencia es en la ordenada en el origen, αi .

¿Por qué ponemos −βi en ver de +βi ?: por analogı́a con regresión logı́stica. En regresión
logistica tenemos una variable W que toma valores 0 ó 1:

p = P r(W = 1|x) logit(p) = β1 x

Si colapsamos las categorı́as de Y (en nuestro caso 1,2,3) en dos:

1, 2 3
|{z}
|{z}
0 1

Entonces:
 
1−p
P r(Y ≤ 2|x) = P r(W = 0|x) = 1 − p ⇒ logit(P r(Y ≤ 2|x) = log = −β1 x
p

Por eso, los βi se interpretan como el cambio en el log-odds de categorı́a más alta frente a
otra más baja, asociado con un incremento en una unidad de Xi , manteniendo las demás
variables constantes. Este cambio es el mismo para cualquiera de las categorı́as ya que, por
ejemplo:

logit(p̃k |Xi = a) = αk − β1 X1 − . . . − βi a + . . . + −βm Xm

logit(p̃k |Xi = a + 1) = α2 − β1 X1 − . . . − βi (a + 1) + . . . + −βm Xm
logit(p̃k |Xi = a + 1) − logit(p̃k |Xi = a) = −βi

9
5 / 64 6 / 64

Proportional odds model

y esto es cierto para cualquier categorı́a k.
Consider a logistic regression model for each logit:
El modelo se llama de odds-proporcionales ya que el odds para cualquier k es proporcio-
nal a logit(θij1 ) = α1 + x0ij β1 None vs. Some/Marked
−β1 X1 +...−βm Xp
=e
0
logit(θij2de
Si una cierta combinación )= α2 + xexplicativas
variables ij β2 None/Some vs.deMarked
duplica el odds estar en la categorı́a
1, también dobla el odds de estar en la categorı́a 2, etc. Los
Proportional odds assumption: regression functions are parallel odds solo se diferencian en el
on the logit
αi
punto de partida
scale i.e., β (cuando
= β . todas las covariables son cero), es decir, en e .
1 2
Proportional Odds Model
Pr(Y>1)
4
1.0 Pr(Y>1)

Pr(Y>2) 3

0.8 Pr(Y>2)
Pr(Y>3) 2
Probability

Pr(Y>3)

Log Odds
1
0.6

0.4
-1

0.2 -2

-3
0.0
-4
0 20 40 60 80 100 0 20 40 60 80 100
Predictor Predictor

7 / 64 8 / 64

Los modelos de regresión ordinal se puede interpretar a partir de variables latentes.

Supongamos que la variable que observamos, Y , proviene de una variable no observada Z,
la cual hemos discretizado:
cj−1 < Z ≤ cj ⇒ Y = j j = 1, . . . , k + 1
La variable Z, depende de covariables:
Z = β1 X1 + . . . + βm Xm + 
donde  tiene una determinada función de distribución F , de modo que:
P r(Y ≤ j) = P r(Z ≤ αj ) = F (αj − β1 X1 + . . . − βm Xm )
y
F −1 (pi + . . . pj ) = αj − β1 X1 + . . . − βm Xm
Si F −1 = logit, tenemos el modelo de odds-proporcionales.

Es importante darse cuenta de que para k = 1 estamos en el caso de una multinomial.

Los OR se calcuları́an:
   
P t(Y ≤ k|X = x1 ) P t(Y ≤ k|X = x0 )
log − = (αk − βx1 ) − (αk − βx0 )
P t(Y > k|X = x1 ) P t(Y > k|X = x0 )
= β11 (x1 − x0 ) ⇒ OR = e−β11 (x1 −x0 )

10
En general un coeficiente positivo indica un incremento en la posibilidad de que un individuo
con valores más altos en la variable explicativa esté en una categorı́a superior de intención
de estudios de postgrado, y un coeficiente negativo indica un descenso la en la posibilidad
de que un individuo con valores más altos en la variable explicativa esté en una categorı́a
superior intención de estudios de postgrado.

2.2. La función polr en R

Esta función forma parte de la librerı́a MASS. Es importante tener en cuenta que esta
función usa una parametrización que implica que los odds ratio son de estar en una categorı́a
superior

Vemos si las variables individuales son significativas en el modelo univariante:

library(MASS)
ord1=polr(apply~pared,postgrado)
summary(ord1)
Coefficients:
Value Std. Error t value
pared1 1.127 0.2634 4.28

Intercepts:
Value Std. Error t value
unlikely|somewhat likely 0.3768 0.1103 3.4152
somewhat likely|very likely 2.4519 0.1826 13.4302
Los parámetros se pueden interpretar de manera similar al caso de regresión logı́stica, pero
ahora hay dos transiciones en vez de una (que serı́a el caso de una variable dicotómica).
El parámetro para pared es positivo lo que indica que el hecho de que los padres tengan
estudios aumenta la posibilidad de hacer estudios de postgrado
ord2=polr(apply~public,postgrado)
summary(ord2)

ord3=polr(apply~gpa,postgrado)
summary(ord3)

ord0=polr(apply~1,postgrado)

anova(ord0,ord1)
anova(ord0,ord2)
anova(ord0,ord3)
El parámetro de gpa también es positivo, lo que indica que la nota media de expediente
predispone a elegir los estudios, y el hecho de que el alumno esté en una Universidad pública
o privada no afecta a la elección de estudios de postgrado. Si introducimos las variables
juntas:

11
ord4=polr(apply~pared+gpa,postgrado)
Antes de seguir interpretando el modelo, hemos de comprobar que se satisface las condicio-
nes del modelo de odds proporcionales, es decir, que si creáramos dos variable dicotómicas
correspodientes a las dos transiciones y a justáramos las variables, los coeficientes en ambos
casos serı́a muy similares. Para comprobarlo podemos utilizar la librerı́a VGAM:
m0 <- vglm(ordered(apply)~pared+gpa,family=cumulative(parallel=T),
data=postgrado)
m1 <- vglm(ordered(apply)~pared+gpa,family=cumulative(parallel=F),
data=postgrado)

test.po <- 2*logLik(m1)-2*logLik(m0)

df.po <- length(coef(m1))-length(coef(m0))
1-pchisq(test.po,df=df.po)
Por lo que podemos aceptar la hipótesis de odds proporcionales.

Calculamos los odds y sus intervalos de confianza:

summary(ord4)
Coefficients:
Value Std. Error t value
pared1 1.0457 0.2656 3.937
gpa 0.6042 0.2539 2.379

Intercepts:
Value Std. Error t value
unlikely|somewhat likely 2.1763 0.7671 2.8370
somewhat likely|very likely 4.2716 0.7922 5.3924

ci=confint(ord4)
exp(cbind(OR=coef(ord4),ci))
OR 2.5 % 97.5 %
pared1 2.845412 1.693056 4.806454
gpa 1.829873 1.115180 3.022320

Entonces, la posibilidad de pasar de nada probable a algo o muy probable, es 2.8 veces mayor
para los alumnos cuyos padres tienen estudios que para los que no los tienen (también serı́a
2.8 veces más probable para de nada o algo probable a muy probable). Algo similar ocurre con
cada unidad que aumenta la nota del expediente. Las ordenadas en el origen corresponderı́a
a los umbrales de corte de la variable latente.

Podemos obtener y dibujar las probabilidades ajustadas:

plot(effect("pared", ord4, style=’stack’))
plot(effect("gpa", ord4, style=’stack’))

12
pr=predict(ord4,postgrado,type="p")

pared effect plot gpa effect plot

very likely very likely
somewhat likely somewhat likely
unlikely unlikely

0.8 0.8
apply (probability)

apply (probability)
0.6 0.6

0.4 0.4

0.2 0.2

0 1 2.0 2.5 3.0 3.5 4.0

pared gpa

No hay métodos equivalentes a los vistos anteriormente para la diagnosis del modelo, con
lo cual tendrı́amos que hacer las regresiones logı́sticas por separado para comprobar si se
cumplen las hipótesis.

2.3. Ejercicios
2.4. Artritis
Los datos corresponden a un estudio sobre el efecto de un determinado tratamiento para
la artritis. La variable respuesta Improved tiene tres categorı́as: None, Some, Marked, y hay
tres variable explicativas: Treatment (Placebo, Tratado), age, sex (Mujer Hombre). Para
obtener los datos:

library(vcd)
data("Arthritis")

¿Qué variables son significativas?, ¿hay interacciones?, ¿cómo interpretas los coeficientes?

3. Regresión de Poisson
Gran cantidad de datos son recogidos por cientı́ficos, médicos, empresas como datos de
conteo. En general, este tipo de datos aparecen en cuantro formatos distintos:

Datos como frecuencias, donde contamos cuantas veces ocurre algo, pero no sabemos
cuantas veces no ocurre.

Datos como proporciones, donde ambos, el número de veces que ocurre algo, y el tamao
toñtal del grupo es conocido.

13
 Datos por categorı́as, donde la varible cuenta cuántos individuos hay en cada nivel de
una variable categórica.

Datos binarios, donde recigemos la presencia o ausencia de una caracterı́astica.

Los datos como proporciones y datos binarios han sido tratados en el capı́tulo de regresión
logı́stica.

Hay cuatro razones por las que es erroneo utlizar un modelo de regresión normal
para datos de conteo :

1. Puede dar lugar a predicciones negativas.

2. La varianza de la variable respuesta no es independiente de la media.

3. Los errores no siguen una distribución Normal.

4. Los ceros que aparecen en la variable respuesta dan problemas a la hora de transformar
la variables.

3.1. La distribución de Poisson

La principal diferencia entre la distribución de Poisson y la Binomial, es que, aunque am-
bas cuentan el número de veces que ocurre algo, en la distribución de Poisson no sabesmos
cuantas veces no ocurrio, y en la Binomial sı́ lo sabemos.

Supongamos que estamos haciendo un estudio sobre cuantas larvas de insectos hay en cier-
tos árboles, los datos de los que disponemos corresponden al número de larvas por hoja (y).
Habrá hojas que no tengan ninguna, y otras que tenga hasta 5 ó 6. Si el número medio de
larvas por hoja es µ, la probabilidad de observar x larvas por hoja viene dada por:
eµ µy
P (y) =
y!
y µ serı́a el número de larvas, por la probabilidad de que una hoja tenga una larva. Implı́ci-
tamente, lo que estamos haciendo es una aproximación:

µ = np

para n grande y p pequeño. Es decir, que una distribución de Poisson se obtiene a partir
de una Binomial con p pequeño y n grande y donde µ = np. Veremos más adelante que
si utilizamos una regresión logı́stica con datos agrupados obtenemos los mismo coeficientes
que si utilizamos una regresión de Poisson, pero la interpretación de ambas sólo es similar
cuando la prevalencia es baja.

El modelo que estamos ajustando es:

log(y) = β0 + β1 X1 + . . . + βK Xk

14
de esta forma nos aseguramos que los valore sajustados para y sean positivos, y eβi representa
un efector multiplicador del i-ésimo predictor sobre la media, ya que al incrementar Xi en
una unidad, la media queda multiplicada por un factor eβi . Una de las ventajas de utilizar
este tipo de modelo es que, en general, en el caso de datos de conteo, los predictores son
multiplicativos y no aditivos, es decir, observamos conteos pequeños para efectos pequeños
y conteos grandes para efectos grandes (o viceversa).

En cuanto a las medidas de bondad de ajuste:

1. Test de la razón de verosimilitud. Sin embargo, hay que tener cuidado con este test de
bondad de ajuste, ya que sólo es válido cuando el número de casos en la mayorı́a de
los estratos es > 5, cosa que a veces no ocurre.
2. Para este tipo de datos utilizamos los residuos estandarizados:
µiobs − µ̂i
ri = √
µ̂i
estos residuos siguen una distribución N (0, 1) por lo tanto la mayorı́a de ellos deben
estar en el intervalo (−2,5, 2,5).

3.2. Regresión de Poisson con variables categóricas y conti-

nuas
Un experimento agrı́cola consta de 90 parcelas de pradera, de 25m2 cada una, que se
diferencian entre sı́, en biomasa, pH del suelo, y la riqueza en especies. Es sabido que
la riqueza de especies decrece cuando aumenta la biomasa, pero la cuestión de interés
aquı́ es si esa disminución cambia con el pH del suelo. Las parcelas se clasificaron de
acuerdo a tres niveles de pH: alto (=2), medio (=1) y bajo (=0). La variable respuesta
es el número de especies por parcela, y las variable predictora son la biomasa media
medida en el mes de Junio, y la variable categórica correspondiente al pH del suelo.

Comenzamos por dibujar los datos:

points(Biomasa[pH==0],Especies[pH==0])
points(Biomasa[pH==1],Especies[pH==1],pch=2,col=2)
points(Biomasa[pH==2],Especies[pH==2],pch=3,col=3)
legend(4,45,legend=c("Bajo","Medio","Alto"),pch=c(1,2,3),col=c(1,2,3))

Viendo el dibujo, podemos tener la tentación de ajustar un modelos de regresión lineal,

pero serı́a erroneo ya que predecirı́a valores negativos para el número de especies si el
valor para la variable Biomasa es alto. Sin embargo, podemos dibujar las lı́neas a modo
de análisis exploratorio:

abline(lm(Especies[pH==0]~Biomasa[pH==0]))
abline(lm(Especies[pH==1]~Biomasa[pH==1]))
abline(lm(Especies[pH==2]~Biomasa[pH==2]))

15
 ● Bajo
Medio

40
Alto

30
Especies

20
● ● ●
●
●
● ●
● ● ●
● ● ●
● ●
● ●
●

10
●
● ●
● ● ●
●

●
● ●
●
●

0 2 4 6 8 10

Biomasa

● Bajo
Medio
40

Alto
30
Especies

20

● ● ●
●
●
● ●
● ● ●
● ● ●
● ●
● ●
●
10

●
● ●
● ● ●
●

●
● ●
●
●

0 2 4 6 8 10

Biomasa

Se puede observar claramente la diferencia en biomasa media para distinto valor del
pH (hay distancia entre las curvas), pero no se aprecia que la pendiente varı́e signifi-
cativamente (las curvas son aproximadamente paralelas). Si hubiéramos ajustado un
modelo de regresión lineal, hubiéramos obtenido que la interacción entre Biomasa y
pH no era significativa. Veamos qué ocurre cuando analizamos los datos utilizando un
modelo de regresión de poisson.

especies1=glm(Especies~Biomasa*pH,family=poisson)
especies2=update(especies1,~.-Biomasa:pH)
anova(especies2,especies1,test="Chi")
Analysis of Deviance Table

Model 1: Especies ~ Biomasa + pH

Model 2: Especies ~ Biomasa * pH
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 86 99.242
2 84 83.201 2 16.040 0.0003288

Lo que nos indica que asumir que descenso en número de especies cuando aumenta la
biomasa es similar sea cual sea el pH, es erroneo.

summary(especies1)
Coefficients:

16
 Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.95255 0.08240 35.833 < 2e-16 ***
Biomasa -0.26216 0.03803 -6.893 5.47e-12 ***
pH1 0.48411 0.10723 4.515 6.34e-06 ***
pH2 0.81557 0.10284 7.931 2.18e-15 ***
Biomasa:pH1 0.12314 0.04270 2.884 0.003927 **
Biomasa:pH2 0.15503 0.04003 3.873 0.000108 ***

El modelo ajustado es el siguiente:

log(Especies) = 2,95−0,261×Biomasa+0,484pHmedio +0,815×pHalto +0,123×pHmedio ×Biomasa+0,1

La interpretación de los parámetros serı́a la siguiente:

Un cambio en una unidad de Biomasa supodrı́a número de especies exp(−0,262) =

0,769 menor en suelos con pH bajo
Un cambio en una unidad de Biomasa supodrı́a número de especies exp(−0,262 +
0,123) = 0,87 menor en suelos con pH medio
Un cambio en una unidad de Biomasa supodrı́a número de especies exp(−0,262 +
0,155) = 0,898 menor en suelos con pH alto

En el siguiente gráfico vemos las curvas ajustadas, y se aprecia cómo el utilizar un

modelo de poisson resuelve el problema de los valores positivos, al ajustar una curva
exponencial en vez de una recta. Si queremos hacer un análisis de residuos:

● Bajo
Medio
40

Alto
30
Especies

20

● ● ●
●
●
● ●
● ● ●
● ● ●
● ●
● ●
●
10

●
● ●
● ● ●
●

●
● ●
●
●

0 2 4 6 8 10

Biomasa

par(mfrow=c(2,2))
plot(especies1)

par(mfrow=c(1,2))
plot(predict(especies1),residuals(especies1))
qqnorm(residuals(especies1))

17
 Normal Q−Q Plot

● ●

3
● ●

2
residuals(especies1)
●● ● ●● ●

Sample Quantiles
● ● ●
● ●● ●●●●
● ● ●●

1
● ● ●● ● ● ●
●
●
●●
●
● ● ● ●●
●
● ●
● ● ●● ● ● ●●
●●
●
●
● ● ●
●
● ● ● ● ● ●
●
●●
●
●● ● ● ● ● ● ●●
●
●●
●
● ● ●●
●

0
● ● ●
●
●●● ● ● ●
●●
●
●
● ● ● ●●
● ●●
● ● ● ● ●● ● ● ● ●
●
●●
●●
●
●
●
● ● ● ● ●●
●
●
● ● ●● ● ●●
●●
● ● ●●

−1

−1
●
● ● ● ● ●●●
●
● ● ●●
● ● ● ● ●
●●●
● ●
●
● ●●
● ●

−2

−2
● ●

2.0 2.5 3.0 3.5 −2 −1 0 1 2

predict(especies1) Theoretical Quantiles

En el caso de datos de Poisson, la media es igual a la varianza, si los datos no cumplen

esta hipótesis los errores estándar de los estimadores de los parámetros serán incorrec-
tos. En estas ocasiones podemos introducir un parámetro de dispersión, φ, de modo
que V ar(Y ) = φµ (si φ = 1 estamos en el caso Poisson). Dicho parámetro se puede
estimar como:
Deviance
φ̂ =
n−p
En el ejemplo anterior, Deaviance = 83,2 y n − p = 84, por lo que el parámetro es
muy próximo a 1. Podemos utilizar en la función glm, la familia quasipoisson, que
automáticamente estima el parámetros de dispersión. dará lugar a los mismos valores
de los parámetros pero cambiarán los errores estándar.

especies2=glm(Especies~Biomasa*pH,family=quasipoisson)
especies2

3.3. Regresión de Poisson para tasas de incidencia

Podemos usar la regresión logı́stica para estimar la prevalencia (proporción), pero no
podemos usarala para estimar la incidencia, ya que en este caso necesitamos utilizar el
tiempo de exposición al riesgo, ya que en la tasas de incidencia medimos el número de
sucesos en la unidad de tiempo.

Supongamos que la variable Y representa el número de eventos que aparecen con

una tasa λ por unidad de tiempo de exposición y que la exposición es E (E se suele
dar en términos de personas por año o 1000×personas por año) y va a determinar las
unidades en las que se expresa la tasa. Dado que nuestro interés es la tasa, λ, y esta
puede depender de distintas covariables, podrı́amos ajustar el modelo:

log(λ) = β0 + β1 X1 + . . . + Xk βk

Pero es Y quien sigue una distribución de Poisson y no λ. Dado que , y ∼ P (µ) y el

número medio de eventos E[Y ] = µ = λ × E:

log(µ) = log(λ) + log(E) = β0 + β1 X1 + . . . + Xk βk + log(E)

18
Al log(E) se le llama offset y es una cantidad fija y no hay que estimar ningún paráme-
tro para ella. La interpretación de los parámetros va unida al concepto de riesgo relativo.
Por ejemplo, supongamos que tenemos una sola covariable X

log(λ) = β0 + Xβ1

Si llamamos λ0 a la tasa de incidencia cuando la covariable toma un valor xo , y λ1 a

la tasa cuando X = x0 + 1, entonces:
 
λ1
log = log(λ1 ) − log(λ0 ) = β1
λ0
λ1
RR = = e β1
λ0
es decir, eβ1 es el riesgo relativo entre la población donde X = x0 y X = x1 , y eβ0 es
la tasa de incidencia cuando todas las variables son cero.

3.4. Ejemplo
Los datos corresponden al registro de incidencias de diabetes mellitus en niños de 0 a
14 años de la Comunidad de Madrid. Las variables son:

Variable Valores
sexo 1=Varón
2=Mujer
edad 1=0-4
2=5-9
3=10-15
mes (numérica del 1 al 12)
periodo 1997-2003
casos número de casos por estrato
poblacio número de personas en el estrato
estacion 1=frı́a (Octubre a Marzo)
2= caliente (Abril a Septiembre)

El objetivo del estudio era:

a) Ver si hay cambios en la incidencia en los 7 años de registro globalmente y por

grupos de edad y sexo
b) Explorar la posibilidad de que la incidencia sea diferente en las diferentes esta-
ciones del año, ası́ como explicar si hay variaciones significativas en la incidencia
según el mes.

diabetes=read.table("diabetes.txt",header=TRUE)
names(diabetes)

19
diabetes$periodo=factor(diabetes$periodo)
diabetes$estacion=factor(diabetes$estacion)
diabetes$sexo=factor(diabetes$sexo)
diabetes$mes=factor(diabetes$mes)
diabetes[1:10,]
edad sexo mes periodo casos poblacion estacion
1 1 2 12 1997 2 110324 1
2 1 2 8 1997 1 110324 2
3 2 2 9 1997 0 120910 2
4 3 2 1 1997 1 146712 1
5 1 1 11 1997 1 116134 1
6 2 2 10 1997 3 120910 2
7 3 1 5 1997 3 154741 2
8 1 2 6 1997 0 110324 2
9 1 2 10 1997 0 110324 2
10 3 2 3 1997 3 146712 1

λ = y/E
por ejemplo para niños, entre 0 y 4 años, en el mes de Diciembre del año 1997, y = 2
y E = 110324, por lo tanto:
2
λ= = 0,0000181 o 1.81 cada 100000 personas
110324
El uso de la función glm() es similar a los ejemplos anteriores, aunque con la peculia-
ridad de que hemos de especificar E, para especificarlo utilizamos offset():

logpobla=log(poblacion)
diabetes1=glm(casos~periodo+offset(logpobla),family=poisson)
summary(diabetes1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -11.21951 0.08944 -125.438 <2e-16 ***
periodo1998 0.02386 0.12649 0.189 0.850
periodo1999 0.04955 0.12599 0.393 0.694
periodo2000 0.13387 0.12391 1.080 0.280
periodo2001 -0.15710 0.13238 -1.187 0.235
periodo2002 0.08106 0.12369 0.655 0.512
periodo2003 0.03106 0.12527 0.248 0.804
---
Null deviance: 686.27 on 503 degrees of freedom
Residual deviance: 680.44 on 497 degrees of freedom

Podemos ver que el año no es una variable significativa. ¿Cómo usarı́as la función
anova() para comprobarlo?.

¿Qué estarı́amos calculando si hacemos:

20
 exp(0.04955-0.02386)

Si queremos obtener los resultados como RR escribimos:

> exp(diabetes1$coeff)
(Intercept) periodo1998 periodo1999 periodo2000 periodo2001 periodo2002
1.340994e-05 1.024145e+00 1.050797e+00 1.143247e+00 8.546228e-01 1.084431e+00
periodo2003
1.031550e+00
> exp(confint(diabetes1))
Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) 1.119437e-05 0.0000159
periodo1998 7.990119e-01 1.3127114
periodo1999 8.206769e-01 1.3456605
periodo2000 8.968513e-01 1.4585479
periodo2001 6.584309e-01 1.1071044
periodo2002 8.511035e-01 1.3829685
periodo2003 8.068865e-01 1.3192888

diabetes1.resid=rstandard(diabetes1)
summary(diabetes1.resid)

La tasa de incidencia ajustada

diabetes1.ajustados=predict(diabetes1,type="response")
tasa1.ajustada=diabetes1.ajustados/poblacion

Ahora estamos listos para responder a las cuestiones que se planteaban al

principio
¡¡MANOS A LA OBRA!!

21