1.- ¿Qué son y por qué hay que utilizar GLM?

Regresión lineal: ¿puede ser la variable de interés predicha por un conjunto de variables explicativas? Entonces, ¿por qué es necesario utilizar otro tipo de modelos?

Para poder utilizar regresión lineal es necesario que la variable respuesta sea continua, y cumpla las hipótesis estándar del modelo lineal (datos Normales, varianza constante, etc.) Si la variable de interés es, por ejemplo binaria e ignoramos este hecho, lo que hacemos es ajustar este modelo:

\[p=Pr(\text{ocurra algo})=\beta_0+\beta_1X_1+\dots+\beta_kX_k\] Al estimar los parámetros usando el procedimiento se podría cometer dos errores graves:

  1. Los valores predichos de la probabilidad podrían estar fuera del intervalo (0,1).
  2. Los procedimientos inferenciales para estudiar qué variables son significativas están basados en la hipótesis de que los datos vienen de una distribución Normal, cosa que no es cierta con datos binarios.

Motivaremos la necesidad de los GLMs usando dos conjuntos de datos donde la distribución de la variable respuesta tiene la propiedad de que la respuesta media y la varianza están relacionadas.

1.1.- Ejemplo de motivación. Percepción de salud

Con datos simulados sobre percepción del estado de salud vamos a ajustar un modelo de regresión lineal simple entre la variable dicotomica que representa el estado de salud y el peso:

set.seed(1234)
peso<-c(rnorm(300,60,8),rnorm(40,100,27))
set.seed(1234)
percepcion<-rbinom(340,1,0.5)
ejemplo=lm(percepcion~peso)
y.ajustado=ejemplo$fitted
plot(peso,percepcion)
abline(ejemplo,col=2)

En 1944, Berkson utilizó por primera vez la regresión logística como una forma de solucionar el problema, ya que la función logit hace que en vez de estar trabajando con valores de la variable respuesta entre (0,1) estemos trabajando con una variable respuesta que puede tomar cualquier valor. No fue hasta 1972 cuando John Nelder introdujo los modelos lineales generalizados (GLM), de ahí que en general se considere la regresión logística como algo distinto a los GLM, cuando lo que ocurre es que tanto la regresión múltiple como la logística, de Poisson, ordinal, etc., son casos particulares de un GLM.

1.2- Ejemplo de motivación. Experimento de turbinas

Se realiza un experimento para determinar la relación entre el tiepo de uso de unas turbinas y el número de fisuras que aparecen. Los datos siguen una distribución de Poisson, por lo tanto, la media es igual a la varianza y la restricción ahora es que los datos ajustados han de ser positivos. Si ajustamos lm(fisuras~horas), obtenemos que el valor ajustado para la primera observación es negativo.

horas<-c(400,1000,1400,1800,2200,2600,3000,3400,3800,4200,4600)
fisuras<-c(0,212,66,511,150,351,378,78,748,840,756)
turbinas<-data.frame(horas,fisuras)
turbinas
##    horas fisuras
## 1    400       0
## 2   1000     212
## 3   1400      66
## 4   1800     511
## 5   2200     150
## 6   2600     351
## 7   3000     378
## 8   3400      78
## 9   3800     748
## 10  4200     840
## 11  4600     756
lm(fisuras~horas)$fitted
##         1         2         3         4         5         6         7         8 
##  -1.47929 101.17751 169.61538 238.05325 306.49112 374.92899 443.36686 511.80473 
##         9        10        11 
## 580.24260 648.68047 717.11834

1.3.- Breve nota sobre las componentes del MLG

Tal y como hemos estudiado, los GLM es establecen una relació lineal no entre la media de la variable respuesta y los predictores (como en la regresión lineal), sino entre una funció de la media de variable respuesta y los predictores, es decir: \[\boldsymbol{\eta}=g(E[Y])=\beta_0+\beta_1X_1+\dots+\beta_rX_r\] según de que tipo sea la varible \(Y\) , así será la función \(g\). Es decir, un GLM tiene 3 componentes:

  1. La variable respuesta \(Y\). Para poder utilizar un GLM, la distribución de \(Y\) ha de pertenecer a la Familia Exponencial, es decir, su función de densidad ha de poder escribirse como: \[f(y;\boldsymbol{\theta},\phi)=\exp\lbrace \frac{y\boldsymbol{\theta}-b(\boldsymbol{\theta})}{a(\phi)}+c(y,\phi) \rbrace\] donde, \(a(\cdot)\), \(b(\cdot)\) y \(c(\cdot)\) son funciones conocidas. El parámetro \(\boldsymbol{\theta}\) se denomina parámetro canónico de localización y \(\phi\) es un parámetro de dispersión.

  2. Las variables predictoras \(X_i, i=1,\dots,k\).

  3. La funcion que relaciona la media, \(E[Y ]\), con las variables predictoras \(X\). En el caso del modelo de regresión ordinaria \(\boldsymbol{\mu}=\boldsymbol{\eta}\), por tanto la función de enlace es la identidad. La función de enlace canónica transforma la media en el parámetro canónico \(\boldsymbol{\theta}.\) Es decir, si \(g(E[Y])=\boldsymbol{\theta},\) entonces \(g\) es función de enlace canónico. La siguiente tabla muestra las funciones de enlace canónico para las distribuciones más comunes:

Funciones de enlace canónico en GLM
Funciones de enlace canónico en GLM

¿Qué diferencia hay entre usar una función de enlace y una transformación?

Concretamente, el modelo de regresión logística que vamos a utilizar es un GLM en la que la variable respuesta es binaria (sigue una distribución de Bernoulli), y la funció de enlace es el logit ya que:

\[g(E[Y]=p)=\log(\frac{p}{1-p})=\beta_0+\beta_1X_1+\dots+\beta_rX_r\]

A continuación ilustraremos el uso en la práctica de la metodología relacionada con los GLM en las clases de teoría, para el modelo de la regresión logística.

2.- Regresión Logística

Vamos a ajustar un modelo de regresión logística a los datos birthwt incluidos en el paquete MASS que, que contiene información sobre el peso de 189 recién nacidos y 10 variables clínicas de sus madres. El objetivo es describir como cambia el bajo peso de los recién nacidos (\(Y\), variable respuesta) al nacer atendiendo la raza, ser o no madre fumadora, el peso de la madre en el último ciclo menstrual o la edad entre otros factores.

library(MASS)
data("birthwt") # Cargamos los datos
dim(birthwt)
## [1] 189  10
head(birthwt)
##    low age lwt race smoke ptl ht ui ftv  bwt
## 85   0  19 182    2     0   0  0  1   0 2523
## 86   0  33 155    3     0   0  0  0   3 2551
## 87   0  20 105    1     1   0  0  0   1 2557
## 88   0  21 108    1     1   0  0  1   2 2594
## 89   0  18 107    1     1   0  0  1   0 2600
## 91   0  21 124    3     0   0  0  0   0 2622
attach(birthwt) # Permite trabajar con los nombres de las variables del data.frame directamente
#-----------------------------------------------------------------
#  low Peso al nacer (0 = Peso $>$ 2500g,
#                     l = Peso $<$ 2500g)
#  age Edad de la madre en a\~nos
#  lwt Peso de la madre en el \'ultimo ciclo menstrual
#  race Raza (1 = Blanca, 2 = Negra, 3 = Otra)
#  smoke Fumadora durante el embarazo (1 = S\'i, 0 = No)
#-----------------------------------------------------------------

Se defina \(y_i=1\) un valor de la variable respuesta cuando el peso del recién nacido es bajo e \(y_i=0,\) en caso contrario. Es decir, la variable \(Y_i\) toma el valor 1 con probabilidad \(p_i\) y o con probabilidad \(1-p_i\) de acuerdo a una distribución de Bernoulli:

\[Pr(Y_i=y_i)=p_i^{y_i}(1-p_i)^{1-y_i},\quad y_i=1,0.\]

En algunas ocasiones, aunque no en este caso, los sujetos estudiados podrían venir dados atendiendo a los factores de interes clasificadosen \(k\) grupos.Por ejemplo, los sujeros podrían clasificarse en 12 grupos atendiendo a ser o no madre fumadora y la raza:

ftable(list(low,smoke,race))
##         x.3  1  2  3
## x.1 x.2             
## 0   0       40 11 35
##     1       33  4  7
## 1   0        4  5 20
##     1       19  6  5

En este caso, se habla de datos agrupados y la formulación cambia ligeramente. Sea \(n_i\) el número de observaciones del grupo \(i\), ahora \(y_i\) denota el número de recién nacidos con bajo peso en el grupo \(i-\)ésimo, donde la variable ahora sigue una distribución binomial de parámetros \(n_i\) y \(p_i\): \[Pr(Y_i=y_i)=\binom{n_i}{y_i}p_i^{y_i}(1-p_i)^{n_i-y_i},\quad y_i=0,1,\dots,n_i.\] La formulación vista en primer lugar y con la que trabajaremos en este ejemplo es un caso particular de la formulación para datos agrupados con \(k=n\) y \(n_i=1.\) Desde el punto de vista práctico, hay que resaltar que si los predictores son categóricos y las observaciones son independientes, ambas maneras de analizar los datos son equivalentes. Una ventaja de trabajar con datos agrupados es que, dependiendo del tama~no de los grupos es posible hacer test sobre la bondad de ajuste del modelo.

Cuando trabajamos/modelizamos probabilidades, el primer problema al que nos enfrentamos es que la probabilidad sólo toma valores entre 0 y 1, y si pretendemos relacionarla directamente con los predictores puede que nos salgamos de dicho intervalo. Una solución sencilla es transformarla y relacionar esta modificació con los predictores. La transformación más común es el o y a este modelo se le conoce como medelo de regresión logística: \[\textit{logit}(p_i)=\ln(\frac{p_i}{1-p_i})=\eta_i=\beta_0+\beta_1X_{1i}+\dots+\beta_rX_{ri}\] de forma que cuando la probabilidad se aproxima a 0 el lo hace a \(-\infty\) y cuando la probabilidad se acerca a 1, el lo hace a \(+\infty\). Si la probabilidad es 0.5, el , ventaja, es 1 y el vale 0. negativos corresponden a probabilidades inferiores a 0.5, y viceversa. Función

La función es biyectiva, a su inverse se le conoce como y permita calcular la probabilidad a partir del : \[p_i=\frac{e^{\eta_i}}{1+e^{\eta_i}}.\] Por tanto, en el modelo de regresión logística, donde el de la probabilidad sigue un modelo lineal, la probabilidad es una función no lineal de los predictores, en concreto: \[p_i=\frac{e^{\beta_0+\beta_1X_{1i}+\dots+\beta_rX_{ri}}}{1+e^{\beta_0+\beta_1X_{1i}+\dots+\beta_rX_{ri}}},\] siendo difícil expresar el efecto que tiene en la probabilidad un cambio en las variables predictoras, esta es una de las razones por las que se trabaja directamente con el .

3.- El modelo de regresión logística como un GLM

El modelo de regresión logística es un GLM, donde la función de enlace canónico es el , es decir, =(),$ \(a(\phi)=1\), \(b(\boldsymbol{\theta})=\ln(1+e^{\boldsymbol{\theta}})\) y \(c(y,\phi)=0\). Por lo que la estimación de los parámetros e inferencia sobre éstos se lleva a cabo a partir de la teoría de verosimilitu vista en clase para los GLMs.

Para estimar los parámetros se maximiza la función de log-verosimilitud (en inglés ). Es ta función representa la verosimilitud de que los datos observados sean de una muestra con una distribución concreta, en nuestro caso Bernoullli: \[\ln(L(\boldsymbol\beta))=\sum\limits_{i=1}^n (y_i\ln(p_i)+(n_i-y_i)\ln(1-p_i))\] La notación \(L(\boldsymbol\beta)\) refleja la dependencia de \(\boldsymbol\theta\) de los parámetros del modelo \(\boldsymbol{\beta}\) (=X), pudiéndose escribir \(p_i\) en función de \(\boldsymbol\beta\). Maximizando (derivando e igualando a cero) esta función, se obtienen los estimadores de los parámetros. Para ello hay que resolver una serie de ecuaciones que no tienen solución anlítica y se utiliza el método iterativo de aproximación de la varosimilitud basado en el algoritmo de Newton-Raphson, ver Apéndice. En concreto, para el ejemplo que nos ocupa, la estimación se calculará mediante el estadístico R.

3.1- Procedimiento para ajustar el modelo de regresión logística en R.

En R la funció que se utiliza para este tipo de modelos es la función glm(). Sus argumentos son: fórmula: similar a la regresión lineal simple (varDep~varInd) familia: binomial, poisson, etc. *link: establece la relación entre la variable respuesta y los predictores: etc.

Por ejemplo si se quiere describir la relación entre el bajo peso al nacer, ser o no madre fumadora y la raza, utilizando esta función:

# Variables categóricas consideradas como tal en R deben definirse como tipo factor
# Se incluyen f-1 variables dummy en la matriz de diseño 
# La categoría de referencia es la primera/baja
smoke<-factor(smoke)
race<-factor(race)

# Procedimiento para ajustar el modelo en R
logistica1=glm(low~smoke+race,family=binomial(link=logit))
summary(logistica1)
## 
## Call:
## glm(formula = low ~ smoke + race, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.8405     0.3529  -5.216 1.83e-07 ***
## smoke1        1.1160     0.3692   3.023  0.00251 ** 
## race2         1.0841     0.4900   2.212  0.02693 *  
## race3         1.1086     0.4003   2.769  0.00562 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 219.97  on 185  degrees of freedom
## AIC: 227.97
## 
## Number of Fisher Scoring iterations: 4

Esta salida incluye la estimación para los parámetros del modelo que comentaremos enseguida. Además, e la última línea de la salida se incluye el número de iteraciones que aparecen son debidas a que es necesario resolver las ecuaciones de verosimilitud en la estimación de los parámetros. El valor del en la primera iteración corresponde a lo que se llama el o modelo nulo, donde todos los parámetros son cero. Si el algoritmo no converge, los estimadores que aparecen no se deben utilizar. El valor se utilizará para comparar modelos anidados. Por último aparece el o criterio de información de Akaike, medida de la calidad relativa de u modelo estadístico que busca la parsimonia entre la bondad de ajuste del modelo y la complejidad del modelo.

En concreto, los valores estimados para los coeficientes del modelo son: \[\ln(\frac{p}{1-p})=-1.84+1.12*\textit{smoke}+1.08*\textit{race2}+1.11*\textit{race3}\] También aparecen los errores estándar asociados con los coeficientes, y son utilizados para contrastar si el parámetro es significativamente distinto de cero y para calcular los intervalos de confianza, que en R se calculan como:

confint(logistica1)
## Waiting for profiling to be done...
##                  2.5 %    97.5 %
## (Intercept) -2.5759256 -1.185417
## smoke1       0.4075730  1.862945
## race2        0.1188796  2.055526
## race3        0.3402930  1.918350

Para poder interpretar los parametros, primero hemos de definir la o ventaja, que es la razón (cociente) entre la probabilidad de experimentar un evento en relación con la probabilidad de no experimentarlo, y que es distinto de la tasa de incidencia. Nuestro modelo está expresado en términos de , y puede escribirse en términos de de la siguiente forma:

exp(coefficients(logistica1))
## (Intercept)      smoke1       race2       race3 
##   0.1587319   3.0526312   2.9567424   3.0300008
exp(confint(logistica1))
## Waiting for profiling to be done...
##                  2.5 %    97.5 %
## (Intercept) 0.07608337 0.3056186
## smoke1      1.50316522 6.4426828
## race2       1.12623429 7.8109464
## race3       1.40535931 6.8097164

Si los datos hubiesen estado agrupados, como se mencionó al inicio de la sección. El ajuste del modelo se realizaría como:

lowBis<-rep(c(0,1),6)
smokeBis<-as.factor(rep(c(0,0,1,1),3))
raceBis<-as.factor(c(rep(1,4),rep(2,4),rep(3,4)))
count<-c(40,4,33,19,11,5,4,6,35,20,7,5)
n<-c(44,44,52,52,16,16,10,10,55,55,12,12)
salud2<-data.frame(cbind(lowBis,smokeBis,raceBis,count,n))
salud2
##    lowBis smokeBis raceBis count  n
## 1       0        1       1    40 44
## 2       1        1       1     4 44
## 3       0        2       1    33 52
## 4       1        2       1    19 52
## 5       0        1       2    11 16
## 6       1        1       2     5 16
## 7       0        2       2     4 10
## 8       1        2       2     6 10
## 9       0        1       3    35 55
## 10      1        1       3    20 55
## 11      0        2       3     7 12
## 12      1        2       3     5 12
salud2=salud2[salud2[,1]==1,]
attach(salud2)
## The following objects are masked _by_ .GlobalEnv:
## 
##     count, lowBis, n, raceBis, smokeBis
summary(glm(cbind(salud2$count,salud2$n-salud2$count)~as.factor(salud2$smokeBis)+
              as.factor(salud2$raceBis),family=binomial(link=logit)))
## 
## Call:
## glm(formula = cbind(salud2$count, salud2$n - salud2$count) ~ 
##     as.factor(salud2$smokeBis) + as.factor(salud2$raceBis), family = binomial(link = logit))
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -1.8405     0.3529  -5.216 1.83e-07 ***
## as.factor(salud2$smokeBis)2   1.1160     0.3692   3.023  0.00251 ** 
## as.factor(salud2$raceBis)2    1.0841     0.4900   2.212  0.02693 *  
## as.factor(salud2$raceBis)3    1.1086     0.4003   2.769  0.00562 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 17.8542  on 5  degrees of freedom
## Residual deviance:  3.1569  on 2  degrees of freedom
## AIC: 31.886
## 
## Number of Fisher Scoring iterations: 4

estuvieran agrupados como aparecen en la tabla aparece al principio del captulo, los comandos en seran:

3.- Interpretación de los parámetros

Los coeficientes estimados para cada variable independiente, representan la tasa de cambio de una función de la variable dependiente (en este caso el ) por unidad de cambio de la variable independiente. Por lo tanto, a la hora de interpretar los parámetros hay que tener en cuenta: tanto la relación funcional entre la variable dependiente y las independientes como la unidad de cambio considerada en las variables independientes. Además la interpretación de los coeficientes dependerá también de qué tipo de variables independientes tengamos: dicotómicas, politómicas o continuas.

3.1.- Variable independiente dicotómica

Comenzaremos con este caso ya que es la base para los demas tipos de variables independientes. Sea \(X\) la variable independiente codificada como 0 o 1 en nuestro modelo de regresión logística, por ejemplo la variable madre fumadora: \[\ln(\frac{p}{1-p})=\beta_0+\beta_1X.\]

La diferencia en el para un individuo con \(X=0\) Y \(x=1\) es \(\beta_1\).

Si en la ecuación anterior despejamos \(p,\) se obtiene: \[p=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}} \text{, y por ende, } 1-p=\frac{1}{1+e^{\beta_0+\beta_1X}}\]

Los distintos valores de las probabilidades para las 4 combinaciones entre la variable dependiente e independiente son:

De acuerdo con esta tabla, la , o ventaja, de que la variable respuesta \(Y\) (bajo peso al nacer) esté presente en los individuos con \(X=1\) (o madre fumadora) se define como: \[\frac{Pr(Y=1|X=1)}{Pr(Y=0|X=1)}=\frac{Pr(Y=1|X=1)}{1-Pr(Y=1|X=1)}\]

De manera análoga, la , o ventaja, de que la variable respuesta \(Y\) (bajo peso al nacer) esté presente en los individuos con \(X=0\) (o madre no fumadora) se define como: \[\frac{Pr(Y=1|X=0)}{Pr(Y=0|X=0)}=\frac{Pr(Y=1|X=0)}{1-Pr(Y=1|X=0)}\]

A partir de estas , podemos definir la , o razón de ventajas, (OR) que es una medida de asociación entre dos variables cualitativas. En términos formales, se define como la probabilidad que una condición de salud o enfermedad (bajo peso al nacer) se presente en un grupo de población expuesto a un agente o factor de riesgo (madre fumadora), sobre la probabilidad de que ocurra en otro grupo sin exposición a dicho riesgo. La OR se define como:\[OR=\frac{Pr(Y=1|X=1)/1-Pr(Y=1|X=1)}{Pr(Y=1|X=0)/1-Pr(Y=1|X=0)}\]

Una OR = 1 implica que no existe asociación entre las variables consideradas, es decir, que la exposición al factor de riesgo no pone a riesgo de enfermar ni protege contra la enfermedad. Los valores mayores de 1 indican que la exposición pone a riesgo de enfermar, mientras que los valores inferiores a 1 indican que ésta protege contra la enfermedad. Por ejemplo, uan OR=3 se interpreta como que los recién nacidos con madres fumadoras tienen 3 veces más posibilidades sufrir bajo peso al nacer.

Utilizando la tabla anterior de probabilidades se tiene: \[OR=\frac{(\frac{e^{\beta_0+\beta1}}{1+e^{\beta_0+\beta_1}})/(\frac{1}{1+e^{\beta_0+\beta_1}})}{(\frac{e^{\beta_0}}{1+e^{\beta_0}})/(\frac{1}{1+e^{\beta_0}})}=\frac{e^{\beta_0+\beta_1}}{e^{\beta_0}}=e^{\beta_1}\]

Nótese que la OR coincide con la exponencial del del coficiente \(\beta_1\) del modelo de regresión logística.

Por ejemplo, supongamos que se quiere estudiar la relación entre madre fumadora y bajo peso al nacer:

table(low,smoke)
##    smoke
## low  0  1
##   0 86 44
##   1 29 30
# Pr(Y=1|X=1)=30/74; Pr(Y=0|X=1)=44/74
# Pr(Y=1|X=0)=29/115; Pr(Y=0|X=0)=86/115
OR<-((30/74)/(44/74))/((29/115)/(86/115)) # OR=(86*30)/(29*44)
OR
## [1] 2.021944
logistica2<-glm(low~smoke,family=binomial(link=logit))
summary(logistica2)
## 
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0871     0.2147  -5.062 4.14e-07 ***
## smoke1        0.7041     0.3196   2.203   0.0276 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 229.80  on 187  degrees of freedom
## AIC: 233.8
## 
## Number of Fisher Scoring iterations: 4
exp(logistica2$coefficients[2])
##   smoke1 
## 2.021944
exp(coefficients(logistica2))
## (Intercept)      smoke1 
##   0.3372093   2.0219436
exp(confint(logistica2))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.2177709 0.5070199
## smoke1      1.0818724 3.8005817

Una OR=2.02 se interpreta como que la ocurrencia de bajo peso en recién nacidos con madres fumadoras es 2.02 veces mayor.

Junto con el estimador de la OR , es importante utilizar los intervalos de confianza (IC) para obtener información adicional acerca del parámetro. Los IC para los coeficientes del modelo se calculan asumiendo que la distribución de los parámetros estimados es aproximadamente normal. Entonces un IC para la OR se calcula tomando exponenciales en el IC para los parámetros:\[exp\lbrace \hat{\beta_1}\pm z_{\alpha/2}\hat{se}(\hat{\beta_1})\rbrace\]

exp(coefficients(logistica2))
## (Intercept)      smoke1 
##   0.3372093   2.0219436
exp(confint(logistica2))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.2177709 0.5070199
## smoke1      1.0818724 3.8005817

Que el IC para el OR sea (1.08,3.80) quiere decir que la ocurrencia de bajo peso al nacer en los recién nacidos de madre fumadora en la población de estudio es entre 1.08 y 3.80 veces más problable (que en los recién nacidos con madre no fumadora).

3.2 Variable independiente politómica

En este caso la variable independiente tiene más de 2 categorías. En nuestro ejemplo, lo ilustraremos con la variable .

logistica3=glm(low~race,family=binomial(link=logit))
summary(logistica3)
## 
## Call:
## glm(formula = low ~ race, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1550     0.2391  -4.830 1.36e-06 ***
## race2         0.8448     0.4634   1.823   0.0683 .  
## race3         0.6362     0.3478   1.829   0.0674 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 229.66  on 186  degrees of freedom
## AIC: 235.66
## 
## Number of Fisher Scoring iterations: 4
exp(coefficients(logistica3))
## (Intercept)       race2       race3 
##   0.3150685   2.3275362   1.8892340
exp(confint(logistica3))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.1929856 0.4950281
## race2       0.9255511 5.7746995
## race3       0.9565420 3.7578709

Para el ajuste de este modelo, R crea tres variables dicotómicas (valores 0/1) que indican si la persona es o no de una determinada raza, a estas variables se las conoce como y se incorporan a la matriz de diseño del GLM. Por defecto, R asume la R que la primera variable, asociada a la primera categoría de la variable (categoría de referencia) toma siempre el valor 0 , y compara el resto de grupos con este. En nuestro caso la categoría de referencia es . La OR para $ es la de las mujeres de dividida por la de las mujeres de , 2.33, esto significa que estimamos que la de bajo peso al nacer para mujeres de es 2.33 veces el de las mujeres con . Los I.C. se obtienen siguiendo el mismo procedimiento que en el caso de las variables dicotómicas.

exp(confint(logistica3))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.1929856 0.4950281
## race2       0.9255511 5.7746995
## race3       0.9565420 3.7578709

3.3 Variable independiente continua

Cuando en un modelo de regresión logística utilizamos una variable predictora continua, la interpretación del parámetro estimado depende de las unidades en que se mide esta variable. En estos casos, el parámetro \(\beta_1\) representa el cambio en el log-odds o cuando la variable independiente (X) cambia en una unidad, o equivalentemente, \(e^{\beta_1}\) es el cambio en la cuando la variable X cambia en una unidad. Si en vez de cambiar en una unidad, cambia en \(c\) unidades, entonces el cambio vendrá dado por \(e^{c\beta_1}.\)

Por ejemplo, si se desea a establecer la relación entre el peso de la madre (en libras) en la última menstruación () y el bajo peso al nacer:

logistica4=glm(low~lwt,family=binomial(link=logit))# lwt como factor?
summary(logistica4)
## 
## Call:
## glm(formula = low ~ lwt, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.99831    0.78529   1.271   0.2036  
## lwt         -0.01406    0.00617  -2.279   0.0227 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 228.69  on 187  degrees of freedom
## AIC: 232.69
## 
## Number of Fisher Scoring iterations: 4
exp(coefficients(logistica4))
## (Intercept)         lwt 
##   2.7137035   0.9860401

Por tanto por cada por cada unidad que aumenta el peso, la de recién nacido con bajo pesp es 0.986 veces menor. Si aumenta 10 unidades, la sería \(0.986^10=0.868\) veces menor. Podemos ver los resultados en un gráfico:

ajustados4=predict(logistica4,type="response")# type=response, devuelve las probabilidades; type=link logit
plot(lwt,ajustados4)

La curva es decreciente ya que la < 1. Hay que tener cuidado al interpretar los resultados para altos ya que están basados en muy pocas observaciones.

3.4 Variables independientes categóricas y continuas

Hasta ahora hemos visto como ajustar modelos de regresión logística con una sola variable independiente, a estos modelos se les llama univariantes. Este tipo de modelos son válidos en pocas ocasiones, ya que en general, una variable independiente está asociada a otras variables independientes. Por eso es necesario considerar un modelo multivariante para poder compreneder mejor lo que ocurre con los datos. El objetivo es ajustar el efecto de cada variable en el modelo por las otras variables independientes. Por lo tanto, en el modelo de regresión logística, cada coeficiente del modelo proporciona una estimación de la log-odds ajustando por las otras variables incluidas en el modelo.

Por ejemplo, consideremos un modelo para explicar la relación entre bajo peso al nacer con madre fumadora (\(X_1\)) y el peso de ésta en su última menstruación (\(X_2\)): \[\ln(\frac{p}{1-p})=\beta_0+\beta_1X_1+\beta_2X_2\]

logistica5=glm(low~smoke,family=binomial(link=logit))
logistica6=glm(low~smoke+lwt,family=binomial(link=logit))
summary(logistica6)
## 
## Call:
## glm(formula = low ~ smoke + lwt, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.62200    0.79592   0.781   0.4345  
## smoke1       0.67667    0.32470   2.084   0.0372 *
## lwt         -0.01332    0.00609  -2.188   0.0287 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 224.34  on 186  degrees of freedom
## AIC: 230.34
## 
## Number of Fisher Scoring iterations: 4

Por tanto, \(e^0.68\) es la OR para un mismo valor de peso, es decir, la OR que esperaríamos encontrar si la el peso medio en madres fumadoras y no fumadoras fuese similar Nótese que en estos casos, la interpretación siempre debe de tener en cuenta la la presencia de otras variables incluidas en el modelo.

Concretamente, se observa que el parámetro \(\beta_1\) es muy similar en ambos modelos (logistica5 y logistica6), esto es debido a que el peso medio de madres fumadoras y no fumadoras es muy similar

summary(logistica5)
## 
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0871     0.2147  -5.062 4.14e-07 ***
## smoke1        0.7041     0.3196   2.203   0.0276 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 229.80  on 187  degrees of freedom
## AIC: 233.8
## 
## Number of Fisher Scoring iterations: 4
mean(birthwt[birthwt[,"smoke"]==0,"lwt"])#kg=lb/2.2046
## [1] 130.8957
mean(birthwt[birthwt[,"smoke"]==1,"lwt"])
## [1] 128.1351

Por tanto ajustar por madre fumadora no cambia la de bajo peso en recién nacido.

ajustados6=predict(logistica6,type="link")
plot(lwt,ajustados6,type="n",main="sin interaccion")
points(lwt[smoke==0],ajustados6[smoke==0],col=2)
points(lwt[smoke==1],ajustados6[smoke==1],col=4)
legend(180,0, col=c(2,4),pch=c(1,1),c("no fumadora","fumadora"))

3.5.- Confusión e interacción

El término confusión se utiliza para describir una covariable que está asociada al mismo tiempo a la variable respuesta y al factor de interés. Cuando estas dos asociaciones están presentes la relación entre la variable respuesta y el factor están confundidos. Una manera de detectar una variable de confusión, siempre que no haya interacción, es ver si los parámetros estimados para el factor de interés cambian ustancialmente al introducir la covariable.

En el caso de bajo peso de recién nacido, para saber si edad () es un variable de confusión con respecto a madre fumadora, vemos si el valor del parámetro asociado a cambia mucho al introducir edad:

logistica7=glm(low~smoke+age,family=binomial(link=logit))
summary(logistica5)
## 
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0871     0.2147  -5.062 4.14e-07 ***
## smoke1        0.7041     0.3196   2.203   0.0276 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 229.80  on 187  degrees of freedom
## AIC: 233.8
## 
## Number of Fisher Scoring iterations: 4
summary(logistica7)
## 
## Call:
## glm(formula = low ~ smoke + age, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.06091    0.75732   0.080   0.9359  
## smoke1       0.69185    0.32181   2.150   0.0316 *
## age         -0.04978    0.03197  -1.557   0.1195  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 227.28  on 186  degrees of freedom
## AIC: 233.28
## 
## Number of Fisher Scoring iterations: 4

En este caso, hay muy poca diferencia, por lo que podemos considerar que no hay confusión.

Por otro lado, tal y como hemos visto, la asociación entre la covariable peso de la madre en la última menstruación ()) y la variable respuesta es la misma para cada nivel del factor (), por tanto se dice que entonces no hay interacción entre la covariable y el factor. Gráficamente, la ausencia de interacción da lugar a un modelo con dos rectas paralelas, una para cada nivel del factor. Cuando la interacción está presente, la asociación entre el factor y la variable respuesta depende del nivel de la covariable, es decir la covariable modifica el efecto del factor. En este caso las rectas no son paralelas, y se dice que la covariable tiene un efecto modificador. Siendo imposible estimar la OR para el factor sin especificar el valor de la covariable para el que se hará dicha comparación.

Veamos en nuestro ejemplo si existe interacción entre el factor madre fumadora y el peso de éstas en su última regla al estudiar la asociación de ambas con la ocurrencia de recién nacido de bajo peso.

logistica8=glm(low~smoke+lwt+smoke:lwt,family=binomial(link=logit))
ajustados8=predict(logistica8,type="link")
summary(logistica8)
## 
## Call:
## glm(formula = low ~ smoke + lwt + smoke:lwt, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  1.93234    1.29209   1.496   0.1348  
## smoke1      -1.51089    1.61737  -0.934   0.3502  
## lwt         -0.02389    0.01039  -2.299   0.0215 *
## smoke1:lwt   0.01757    0.01279   1.373   0.1697  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 222.37  on 185  degrees of freedom
## AIC: 230.37
## 
## Number of Fisher Scoring iterations: 4
plot(lwt,ajustados6,type="n",main="sin y con interaccion")
points(lwt[smoke==0],ajustados6[smoke==0],col=2)
points(lwt[smoke==1],ajustados6[smoke==1],col=4)
points(lwt[smoke==1],ajustados8[smoke==1],col=3,pch=2)
legend(180,0, col=c(2,4,3),pch=c(1,1,2),c("no fumadora","fumadora", "fumadora-interac"))

El gráfico muestra que introducir un término de interacción modifica la pendiente de la recta, indicio de que que hay interacción, es un efecto modificador. ¿Será una variable de confusión?

summary(logistica5)
## 
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0871     0.2147  -5.062 4.14e-07 ***
## smoke1        0.7041     0.3196   2.203   0.0276 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 229.80  on 187  degrees of freedom
## AIC: 233.8
## 
## Number of Fisher Scoring iterations: 4
summary(logistica6) # No considerar la interaccion porque esta modifica la estimacion de los parametros
## 
## Call:
## glm(formula = low ~ smoke + lwt, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.62200    0.79592   0.781   0.4345  
## smoke1       0.67667    0.32470   2.084   0.0372 *
## lwt         -0.01332    0.00609  -2.188   0.0287 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 224.34  on 186  degrees of freedom
## AIC: 230.34
## 
## Number of Fisher Scoring iterations: 4

La estimación del parámetro pasa de 0.67 a 0.7 por lo que no parece que el peso sea una variable de confusión (un decremento/incremento de más de 25% sería un indicio de variable de confusión). En cualquier caso, ambas conclusiones hay que comprobarlas mediante tests de hipótesis que comparen los modelos.

Veamos un último ejemplo de un modelo donde la interacción no es significativa. Como vimos la edad no es una variable de confusión, ¿existe interacción entre el factor madre fumadora y la edad para nuestra variable repuesta?

#logistica7=glm(low~smoke+age,family=binomial(link=logit))
summary(logistica5)
## 
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0871     0.2147  -5.062 4.14e-07 ***
## smoke1        0.7041     0.3196   2.203   0.0276 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 229.80  on 187  degrees of freedom
## AIC: 233.8
## 
## Number of Fisher Scoring iterations: 4
summary(logistica7)
## 
## Call:
## glm(formula = low ~ smoke + age, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.06091    0.75732   0.080   0.9359  
## smoke1       0.69185    0.32181   2.150   0.0316 *
## age         -0.04978    0.03197  -1.557   0.1195  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 227.28  on 186  degrees of freedom
## AIC: 233.28
## 
## Number of Fisher Scoring iterations: 4
logistica9=glm(low~smoke+age+smoke:age,family=binomial(link=logit))
summary(logistica9)
## 
## Call:
## glm(formula = low ~ smoke + age + smoke:age, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.80499    1.02523   0.785   0.4323  
## smoke1      -0.96342    1.51147  -0.637   0.5239  
## age         -0.08288    0.04499  -1.842   0.0655 .
## smoke1:age   0.07308    0.06534   1.118   0.2634  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 226.01  on 185  degrees of freedom
## AIC: 234.01
## 
## Number of Fisher Scoring iterations: 4
ajustados9=predict(logistica9,type="link")
plot(age,ajustados9,type="n",main="con interaccion")
points(age[smoke==0],ajustados9[smoke==0],col=2)
points(age[smoke==1],ajustados9[smoke==1],col=4)
legend(20,-2.5, col=c(2,4),pch=c(1,1),c("no fumadora","fumadora"))

En el gráfico vemos que las rectas a las que da lugar el modelo con interacción no son paralelas, pero más importante el coeficiente asociado a la interacción no es significativo y no debe considerarse. En la práctica, cualquier cambio significativo en el coeficiente del factor de riesgo, sugiere que la variable es de confusión. Si esto ocurre y la interacción no es estadísticamente significativa, la variable ha de salir del modelo. Por otra parte, una variable es un efecto modificador, solo cuando el término de interacción añadido al modelo es estadísticamente significativo.

3.6.- Interpretación de la OR en presencia de interacción

Cuando existe interacción entre un factor de riesgo y otra variable, el parámetro estimado para el factor de riesgo depende del valor de la variable que interactúa con este. Luego, no es posible obtener la OR simplemente exponenciando el valor del parámetro, los pasos a seguir son:

  • Escribir la ecuación del para los dos niveles del factor de riesgo.
  • Calcular la diferencia entre los .
  • Exponenciar el valor obtenido.

Por ejemplo, en un modelo que contiene solo dos variables y su interacción. Llamamos \(F\) al factor, \(X\) a la covariable y \(F\times X\) a la interacción. El para \(F=f\) y \(X=x\) es: \[\ln(\frac{p(f,x)}{1-p(f,x)})=\beta_0+\beta_1f+\beta_2x+\beta_3f\times x,\] si queremos la OR para los niveles de F \(f_1\) \(f_0\): \[\begin{align*} \ln(\frac{p(f_1,x)}{1-p(f_1,x)})&=\beta_0+\beta_1f_1+\beta_2x+\beta_3f_1\times x\\ \ln(\frac{p(f_0,x)}{1-p(f_0,x)})&=\beta_0+\beta_1f_0+\beta_2x+\beta_3f_0\times x\\ \ln(\frac{p(f_1,x)}{1-p(f_1,x)})-\ln(\frac{p(f_0,x)}{1-p(f_0,x)})&=\beta_1(f_1-f_0)+\beta_3 x(f_1-f_0)\\ OR&=exp\lbrace \beta_1(f_1-f_0)+\beta_3 x(f_1-f_0) \rbrace \end{align*}\]

En nuestro modelo con interacción (que asumiremos significativa) entre el peso de la madre en su última menstruación y ser madre fumadora o no, la ventaja de recién nacido de bajo peso de madre fumadora frente a no fumadora para una madre de peso en la última menstruación \(x\) se calcula como:

summary(logistica8) # Asumiremos que sola interaccion es significativa para ilustrarlo
## 
## Call:
## glm(formula = low ~ smoke + lwt + smoke:lwt, family = binomial(link = logit))
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  1.93234    1.29209   1.496   0.1348  
## smoke1      -1.51089    1.61737  -0.934   0.3502  
## lwt         -0.02389    0.01039  -2.299   0.0215 *
## smoke1:lwt   0.01757    0.01279   1.373   0.1697  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 222.37  on 185  degrees of freedom
## AIC: 230.37
## 
## Number of Fisher Scoring iterations: 4
# OR se calcula:
#exp(-1.51089+0.01757*x)#x es un valor del peso de la madre

3.7.- Interpretación de los valores ajustados

En general, en regresión logística, los parámetros y las OR son el centro de interées, aunque en ocasiones los valores ajustados pueden ser también importantes. Calcular los IC para el es sencillo, ya que R da los errores estándar. Para calcular el IC para las probabilidades sÓlo hay que utilizar la relaciÓn entre el y la probabilidad: \[\frac{e^{IC logit}}{1+e^{IC logit}}\] Por ejemplo, supongamos que nos interesa el efecto del peso de la mujer en la última menstruación sobre la ocurrencia de recién nacidos de bajo peso en madres fumadoras.

logistica6=glm(low~smoke+lwt,family=binomial(link=logit))
ajustados6=predict(logistica6, se.fit=TRUE)
names(ajustados6)
## [1] "fit"            "se.fit"         "residual.scale"
L.inf=with(ajustados6,exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit))) # Evalua un expresión de R con la función pasada como argumento
L.sup=with(ajustados6,exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)))
with(ajustados6,plot(lwt[smoke==1],(exp(fit)/(1+exp(fit)))[smoke==1],
                     ylim=c(0,1),ylab="probabilidad"))
points(lwt[smoke==1],L.inf[smoke==1],col=2)
points(lwt[smoke==1],L.sup[smoke==1],col=2)
legend(150,0.9, c("probabilidad","I.C."),col=c(1,2),pch=c(1,1))
abline(v=mean(lwt[smoke==1]),lty=2)
abline(v=135)

# ¿Y en las no fumadoras?

El IC es más ancho en los extremos ya que hay menos observaciones y por lo tanto las estimaciones son menos precisas. La línea vertical corresponde al peso medio de las mujeres fumadoras, donde el IC es más estrecho. Por ejemplo, para un peso de 130 libras, el valor estimado es 0.45 y el IC [0.3,0.55].

4.- Bondad de ajuste

Un aspecto importante en el ajuste de cualquier modelo, en particular en los GLM, es determinar si describe adecuadamente o no los datos observados, para ello utilizaremos el test de la razón de verosimilitudes y el test de Wald

4.1- Deviance. Test razón de verosimilitudes

En el caso de un GLM y en particular en la regresión logística la filosofía detrás de un contraste de bondad de ajuste es commparar los valores observados y predichos por el modelo con y sin la variable explicativa. En un GLM la comparación se hace utilizando la función de verosimilitud. Si definimos un modelo saturado como aquel en que hay tantos parámetros como observaciones, la comparación se plantea a través del estadístico Deviance (D) definido como:

\[D=-2\ln(\frac{\text{verosimilitud del modelo ajustado}}{\text{verosimilitud del modelo saturado}}),\] a la cantidad de dentro de los corchetes se llama razón de verosimilitud. Dicho estadístico mide como de grande es la desviación del modelo ajustado respecto de los datos (modelo saturado). Si
D es grande, el modelo ajustado proporciona un ajuste pobre. Un valor pequeño de
D es indicador de un buen ajuste. Para decidir qué se considera D grande" y quéD pequeño”, es preciso utilizar la distribución en el muestreo del estadístico.

En concreto, cuando el modelo propuesto proporciona un buen ajuste de los datos, entonces \(D\sim\chi^2_n-r}\) asintóticamente, donde \(n\) es el número de datos a ajustar y \(r\) el número de parámetros del modelo. Conocida la distribución es posible realizar contrastes de hipótesis. A este test se le conoce como test de la razón de verosimilitudes.

En el caso particular de que en modelo de regresión logística haya una única variable dicotómica, el modelo saturado es 1 y por lo panto la Deviance viene dada por: \[D=-2\ln(\text{verosimilitud del modelo ajustado})\]

logistica2
## 
## Call:  glm(formula = low ~ smoke, family = binomial(link = logit))
## 
## Coefficients:
## (Intercept)       smoke1  
##     -1.0871       0.7041  
## 
## Degrees of Freedom: 188 Total (i.e. Null);  187 Residual
## Null Deviance:       234.7 
## Residual Deviance: 229.8     AIC: 233.8
anova(logistica2,test="Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: low
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                    188     234.67           
## smoke  1   4.8674       187     229.81  0.02737 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Deviance
1-pchisq(logistica2$deviance,logistica2$df.residual)
## [1] 0.0179593

4.2- Test de Wald

Para determinar la significancia individual de cada uno de los predictores introducidos en un modelo de regresión logística se emplea el estadístico Z y el test Wald chi-test. Este test compara el estimador máximo verosímil del parámetro con su error estándar, bajo la hipótesis nula \(H_0: \beta=0\). La razón entre estas dos cantidades sigue una distribución \(N(0,1)\), pero este test lleva a no rechazar \(h_0\) cuando el parámetro es significativo.

logistica10<-glm(low~lwt+smoke+race,family=binomial)
summary(logistica10)
## 
## Call:
## glm(formula = low ~ lwt + smoke + race, family = binomial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.10922    0.88211  -0.124  0.90146   
## lwt         -0.01326    0.00631  -2.101  0.03562 * 
## smoke1       1.06001    0.37832   2.802  0.00508 **
## race2        1.29009    0.51087   2.525  0.01156 * 
## race3        0.97052    0.41224   2.354  0.01856 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 215.01  on 184  degrees of freedom
## AIC: 225.01
## 
## Number of Fisher Scoring iterations: 4

4.3- \(R^2\) y otras medidas predictivas

Cuando hablamos de regresión lineal el coeficiente de correlación y el de determinación son medidas útiles para saber cómo de bien se ajusta el modelo a los datos. En regresión logística podemos calcular una medida análoga al \(R^2\) definida como:

\[R^2_L=\frac{\text{log-verosimilitud del modelo ajustado}-\text{log-verosimilitud del modelo nulo}}{\text{log-verosimilitud del modelo saturado}-\text{log-verosimilitud del modelo nulo}}\]

donde el modelo nulo es aquel ajustado únicamente con el término del intercept. \(R^2_L\) varía entre 0 y 1. Es igual a 0 cuando el modelo no mejora en nada sobre el modelo nulo, y es igual a 1 cuando el modelo ajusta tan bien como el modelo saturado. Una debilidad es que la escala para la verosimilitud logarítmica puede no ser tan fácil de interpretar como la escala para la variable de respuesta en sí. La medida es principalmente útil para comparar modelos. Siguiendo esta línea pueen calcularse otros estadísticos similares.

library(DescTools)
PseudoR2(logistica1,which="all")
##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##      0.06262906      0.02853892      0.07481672      0.10521367      0.07215258 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##      0.13026279      0.06820674      0.11117739      0.07223878    227.97471054 
##             BIC          logLik         logLik0              G2 
##    240.94169860   -109.98735527   -117.33599810     14.69728565

5.- Comparación y selección de modelos

Aunque un modelo con más parámetros proporcione un ajuste mejor que otro con menos parámetros, cabe plantearse si en realidad todos los parámetros estimados son necesarios. Esta cuestión se puede resolver con tests basados en la deviance y en el estadístico de Wald, siempre y cuando se trate de modelos anidados. Otros estadísticos como el AIC, BIC, permiten la comparación más general de modelos no necesariamente anidados. Tmbién pueden ser interesantes medidas de la capacidad predictiva del modelo como el \(R^2\) o análogos.

5.1- Razón de Verosimilitud. Diferencia de Deviances.

Dos modelos en competencia pueden ser comparados mediante la deviance cuando tienen la misma distribución y función link y sólo difieren en el número de parámetros, es decir, cuando se trata de modelos anidados. Supongamos dos modelos M1 con \(p_1\) parámetros y \(M2\) con \(p_2\) parámetros, de forma que el primero está anidado en el segundo, es decir, \(p_1<p_2\). En esta situación el estadístico para valorar si ambos modelos pueden considerarse iguales es: \[G=-2\ln\frac{\text{verosimilitud M1 }}{\text{verosimilitud M2}}=D(\text{verosimilitud M1 })-D(\text{verosimilitud M2 })\] El estadístico \(G\) se corresponde con la diferencia de Deviance y asintóticamente sigue \(D\sim\chi^2_{p_2-p_1}}.\) asintóticamente

anova(glm(low~smoke,family=binomial),
      glm(low~smoke+race,family=binomial),test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: low ~ smoke
## Model 2: low ~ smoke + race
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       187     229.81                        
## 2       185     219.97  2   9.8299 0.007336 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(glm(low~smoke,family=binomial),
      glm(low~smoke+race+lwt,family=binomial),test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: low ~ smoke
## Model 2: low ~ smoke + race + lwt
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       187     229.81                        
## 2       184     215.01  3    14.79 0.002005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.2- Criterios de información

Como en la regresión lineal podemos usar el criterio de información de Aike (AIC) y el criterio de información de Bayes (BIC) para juzgar el ajuste del modelo. Estos criterios proporcionan una medida del ajuste de un modelo que penaliza al modelo que contiene más variables predictoras, pudiendo así comparar dos modelos.

\[\begin{align*} AIC&=−2\text{log-verosimilitud}+2r\\ BIC&=−2{log-verosimilitud}+2r×\log(n)$ \end{align*}\]

En la comparación de modelos, en ambos casos seleccionaremos el modelo con menor valor en dichos estadísticos. Se debe recordar siempre que el mejor modelo es aquel con mayor log-verosimilitud y menor AIC/BIC.

logLik(logistica10) # Valor de la log-verosimilitud
## 'log Lik.' -107.5073 (df=5)
AIC(glm(low~smoke,family=binomial))
## [1] 233.8046
AIC(glm(low~smoke+race,family=binomial))
## [1] 227.9747
AIC(glm(low~smoke+race+lwt,family=binomial))
## [1] 225.0147
BIC(glm(low~smoke,family=binomial))
## [1] 240.2881
BIC(glm(low~smoke+race,family=binomial))
## [1] 240.9417
BIC(glm(low~smoke+race+lwt,family=binomial))
## [1] 241.2234

5.3- Métodos paso a paso o stepwise

Si usamos el método hacia delante () se comienza con un modelo que incluye solo la constante y entonces añade predictores individuales al modelo basándose en la variable que mejore el AIC o BIC. El ordenador continua mientras ninguno de los predictores restantes haya disminuido el criterio.

El método hacia atrás () usa el mismo criterio pero empieza el modelo incluyendo todas las variables predictoras. R comprueba si alguno de esos predictores puede ser retirado del modelo sin incrementar el criterio de información. Si se puede, esa variable se saca del modelo, y se analizan de nuevo el resto de variables.

Mejor que estos métodos es el de ambas direccioes que empieza como el método hacia delante, con solo la constante, pero cada vez que añadimos una variable, el ordenador comprueba si merece la pena eliminarla.

library(RcmdrMisc)
## Loading required package: car
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
## 
##     Recode
## Loading required package: sandwich
modelo <- stepwise(logistica10, direction='backward', criterion='AIC')
## 
## Direction:  backward
## Criterion:  AIC 
## 
## Start:  AIC=225.01
## low ~ lwt + smoke + race
## 
##         Df Deviance    AIC
## <none>       215.01 225.01
## - lwt    1   219.97 227.97
## - race   2   224.34 230.34
## - smoke  1   223.26 231.26
modelo <- stepwise(logistica10, direction='forward', criterion='AIC')
## 
## Direction:  forward
## Criterion:  AIC 
## 
## Start:  AIC=236.67
## low ~ 1
## 
##         Df Deviance    AIC
## + lwt    1   228.69 232.69
## + smoke  1   229.81 233.81
## + race   2   229.66 235.66
## <none>       234.67 236.67
## 
## Step:  AIC=232.69
## low ~ lwt
## 
##         Df Deviance    AIC
## + smoke  1   224.34 230.34
## + race   2   223.26 231.26
## <none>       228.69 232.69
## 
## Step:  AIC=230.34
## low ~ lwt + smoke
## 
##        Df Deviance    AIC
## + race  2   215.01 225.01
## <none>      224.34 230.34
## 
## Step:  AIC=225.01
## low ~ lwt + smoke + race
modelo <- stepwise(logistica10, direction='forward/backward', criterion='AIC') # No tiene porque llegarse al mismo resultado
## 
## Direction:  forward/backward
## Criterion:  AIC 
## 
## Start:  AIC=236.67
## low ~ 1
## 
##         Df Deviance    AIC
## + lwt    1   228.69 232.69
## + smoke  1   229.81 233.81
## + race   2   229.66 235.66
## <none>       234.67 236.67
## 
## Step:  AIC=232.69
## low ~ lwt
## 
##         Df Deviance    AIC
## + smoke  1   224.34 230.34
## + race   2   223.26 231.26
## <none>       228.69 232.69
## - lwt    1   234.67 236.67
## 
## Step:  AIC=230.34
## low ~ lwt + smoke
## 
##         Df Deviance    AIC
## + race   2   215.01 225.01
## <none>       224.34 230.34
## - smoke  1   228.69 232.69
## - lwt    1   229.81 233.81
## 
## Step:  AIC=225.01
## low ~ lwt + smoke + race
## 
##         Df Deviance    AIC
## <none>       215.01 225.01
## - lwt    1   219.97 227.97
## - race   2   224.34 230.34
## - smoke  1   223.26 231.26

6.- Diagnóstico del modelo

Los residuos se pueden utilizar para explorar la adecuación del ajuste de un modelo respecto a la elección de la función de varianza, la función link y los términos a incluir en el predictor lineal. Los residuos pueden también ser útiles para detectar observaciones influyentes o valores anómalos que requieran una investigación más intensa. Para GLM, requerimos una definición algo más amplia de residuos, que sea aplicable a todas las distribuciones que pueden reemplazar a la Normal, tal y como se menciona en el libro de Agresti..

La especificación incorrecta de un modelo se puede deber a una serie de factores:

elección incorrecta de la distribución de probabilidad, especificación incorrecta de la forma en que la respuesta media cambia con las variables explicativas, ya sea porque: la componente sistemática está mal especificada, o la función link no es apropiada, variables explicativas no incorporadas en el modelo, funciones incorrectas de las variables explicativas en el modelo, incluyendo interacciones desconocidas entre ellas, o no disponibilidad de suficientes funciones diferentes, *dependencia entre las observaciones, por ejemplo, a lo largo del tiempo (autocorrelación). Se definen en este caso los residuos deviance estandarizados que son los utilizados para diagnosticar el modelo, obtenidos como los residuos asociados con el predictor lineal. Sin embargo, dada la estructura de los GLM podemos definir también los residuos de la respuesta que son los obtenidos entre el valor de la respuesta y el proporcionado al deshacer el cambio involucrado en la función link.

# Residuos Chi-cuadrado
residuos.p<-resid(logistica2,"pearson") 
X2<-sum(residuos.p*residuos.p) 
X2
## [1] 189
1-pchisq(X2,logistica2$df.residual)
## [1] 0.4453171
# No sirve variables cuantitativas


# Diágnosis del modelo 

A<-influence.measures(logistica10)$infmat
dr=residuals(logistica10,"deviance")
pr=predict(logistica10,type="response")
plot(pr,dr)

plot(logistica10)

plot(dr,A[,6])

plot(dr,A[,1])

Como se puede observar en el gráfico de ajustados versus residuos la interpretación en este tipo de modelos se hace bastante complicada. Esto es debido a que la variable respuesta sólo toma valores 0 y 1, lo que motiva que el gráfico tenga esa pinta tan extraña, se identifican cuatro observaciones atípicas. Si valoraros las distancias de Cook obtenidas podemos ver que no hay ninguna superior a 1, aunque una observación podría considerarse de influencia.

7.- Predicciones y clasificación con el modelo de regresión logística

Una de las principales aplicaciones de un modelo de regresión logística es que puede ser empleada como una herramienta de clasificación. La siguiente tabla resulta de cruzar la variable respuesta con una variable dicotómica cuyos valores resultan de las probabilidades estimadas en el modelo. Para definir esta variable hemos de decidir un punto de corte c y comparar las probabilidades estimadas con c, si son mayores, la variable vale 1 y si son menores 0. Una posibilidad estándar es considerar c=0.5 , aunque podría ser cualquier otro valor.

En ningún caso, las tablas no deben usarse para comparar modelos ya que depende en gran parte de la distribución de las probabilidades en la muestra en la que están basadas, el mismo modelo evaluado en dos muestras diferentes puede dar lugar a clasificaciones distintas. Si el objetivo del análisis no es la clasificación de sujetos, entonces las tablas de clasificación son sólo el complemento de otras medidas de ajuste.

Una de las medidas más utilizadas para valorar la calidad de la predicción el el área bajo la curva ROC () para ello debemos definir primeramente algunos conceptos: Sensitividad: proporción de verdaderos 1s estimados como 1s (la probabilidad de predecir correctamente un 1): Ss = a/(a + c). Especificidad es la proporción de verdaderos 0s estimados como 0s (la probabilidad de predecir correctamente un 0): Sp = d/(b + d). Tasa de falsos positivos: proporción de verdaderos 0s estimados como 1s (la probabilidad de predecir incorrectamente un 0): F+ = b/(b + d). Tasa de falsos negativos es la proporción de verdaderos 1s estimados como 0s (la probabilidad de predecir incorrectamente un 1): F- = c/(a + c).

Cada vez que elegimos un punto de corte la sensitividad y especificidad cambian, si nuestro objetivo es buscar el punto de corte óptimo desde el punto de vista de la clasificación, entonces buscaremos aquel que maximiza la sensitividad y especificidad. La curva ROC es un gráfico de sensibilidad frente a 1-especificidad (tasa de falsos positivos) para todos los puntos de corte. El área por debajo de esa curva puede tomar valores entre 0 y 1, y nos proporciona una medida de la habilidad del modelo para discriminar entre aquellos individuos que presentan la respuesta de interes frente a los que no la presentan.

En general: Si ROC <0.5 el modelo no ayuda a discriminar Si 0;6 <ROC<0.8 el modelo discrimina de forma adecuada Si 0;8 <ROC <0.9 el modelo discrimina de forma excelente Si ROC <0.9 el modelo discrimina de forma excepcional

pred10<-predict(logistica4,type="response")

table(low,ifelse(pred10<0.1,0,1))
##    
## low   0   1
##   0   4 126
##   1   0  59
table(low,ifelse(pred10<0.2,0,1))
##    
## low   0   1
##   0  17 113
##   1   4  55
table(low,ifelse(pred10<0.3,0,1))
##    
## low  0  1
##   0 51 79
##   1 14 45
table(low,ifelse(pred10<0.339,0,1))
##    
## low  0  1
##   0 87 43
##   1 30 29
table(low,ifelse(pred10<0.4,0,1))
##    
## low   0   1
##   0 120  10
##   1  50   9
library(Epi)
ROC(form=low~lwt+smoke+race,plot="ROC")

8.- Ejercicio propuesto

El transatlántico británico RMS Titanic se hundió después de chocar con un iceberg en su viaje inaugural. La tragedia se saldó con la muerte de 1502 de los 2224 pasajeros y tripulantes. En el siguiente enlace https://web.stanford.edu/class/archive/cs/cs109/cs109.1166/stuff/titanic.csv puedes descargar los datos relativos a los pasajeros, como por ejemplo si sobrevivieron, la clase en la que viajaban, el número de hijos a bordo etc. La descripción de las variables aparece a continuación:

survival - Survival (0 = No; 1 = Yes) class - Passenger Class (1 = 1st; 2 = 2nd; 3 = 3rd) name - Name sex - Sex age - Age sibsp - Number of Siblings/Spouses Aboard parch - Number of Parents/Children Aboard ticket - Ticket Number fare - Passenger Fare cabin - Cabin embarked - Port of Embarkation (C = Cherbourg; Q = Queenstown; S = Southampton) boat - Lifeboat (if survived) *body - Body number (if did not survive and body was recovered)

Si bien hubo algún elemento de suerte involucrado en la supervivencia, parece que algunos grupos de personas tenían más probabilidades de sobrevivir que otros. Analice varios modelos de regressión logística y seleccione uno de ellos para modelizar la supervivencia de los pasajeros. Para el modelo seleccionado incluya la bondad del ajuste, diagnosis del modelo y las principales conclusiones que se puedan extraer. ## Apéndice. Estimación en los modelos GLM