Referencias

Las notas se han elaborado a partir de vario material:

Base teórica

Ya conocemos los modelos de regresión lineal mediante los cuales se explica una variable teórica \(Y\) en función de \(k\) variables \(X_1,...,X_n\) independientes o covariables, suponiendo una relación lineal del tipo \[ E[Y]=\beta_0+\beta_1x_1+...+\beta_k x_k \] donde los \(x_i\) son los valores dados de lax variables \(X_i\). Hemos visto que en este tipo de análisis es que la variable de respuesta \(Y\) es de tipo contínuo con distribución normal.

En la Regresión Logística la variable de respuesta es de tipo dicotómico, tomando pués dos valores, el de éxito y el de fracaso, siendo \(p\) la probabilidad que \(Y\) tome el valor de éxito.

En este tipo de situación, el modelo que describe la relación entre la variable de respuesta \(Y\) con las regresoras, \(X_1,...,X_k\) es \[ \text{ln}\frac{p}{1-p}=\beta_0+\beta_1x_1+...+\beta_kx_k \]

Donde las covariables \(X_j\) pueden ser de tipo cualitativo o cuantitativo.

Despejando la \(p\), la relación anterior quedaría como1

\[ p=\frac{\text{exp}(\beta_0+\sum_{j=1}^k \beta_jx_j)}{1+\text{exp}(\beta_0+\sum_{j=1}^k \beta_jx_j)}= \frac{e^{b_0+\sum_{j=1}^k\beta_jx_j}}{1+e^{b_0+\sum_{j=1}^k\beta_jx_j}} \]

que definimos como Regresión logística, Mientras \(\frac{p}{1-p}\) es la Odds Ratio (OR) y \(\text{log}\frac{p}{1-p}\) se denomina \(\text{logit}\)

Incorporación de las variables cualitativas al modelo

Las variables cualitativas del modelo se tienen que convertir en tantas variables dummy cuantas clases presente la variable menos 1. De esta forma, todas las varibles del modelo logísitco acabarán siendo cuantitativas.

De esta forma, los modelos de regresión logística suelen presentar muchas covariables cualitativas que deben ser convertidas en variables indicadoras (dummy), previamente al análisis.

Ejemplo 1 con R

Datos

library(dplyr)
inf <- read.table("datosr/infarto.txt", header = T)

inf <- inf %>% mutate(across(2:5, ~as.factor(.))) # El dataframe es numérico. Se convierte en "factor" para poder ser procesado por la función "model.matrix"

Trasformación de variables categóricas en variables indicadores (dummy variables)

Con R, podemos llevar a cabo las operaciones de tranformación de varibles categóricas en dummy varibles, mediante la función \(\tt model.matrix()\). En este caso las variables son numéricas porque los valores categóricos se han tranformado en numéricos. P.ej: \(Actividad\)

\[ \begin{align} & \texttt{Actividad1=(1,0,0)} = \text{Nula}\\ & \texttt{Actividad2=(0,1,0)} = \text{Mínima}\\ &\texttt{Actividad3=(0,0,1)} = \text{Moderada}\\ &\texttt{Actividad4=(0,0,0)} = \text{Más que moderada} \end{align} \]

Por lo que acabamos de explicar, \(\tt Actividad4\) no va a salir explícitamente en el modelo.

inf.dummy <- model.matrix(~ Actividad+Infarto+Angina+Presion, data = inf)[,2:7]
head(inf.dummy)
##   Actividad1 Actividad2 Actividad3 Infarto1 Angina1 Presion1
## 1          1          0          0        1       1        1
## 2          0          1          0        1       1        1
## 3          0          1          0        1       1        1
## 4          0          1          0        1       1        1
## 5          1          0          0        1       1        1
## 6          0          0          0        1       0        1

De seguida volvemos a incluir la variable numérica junto con las dummy variables.

valores <- bind_cols(select(inf,1), as.data.frame(inf.dummy))
head(valores)
##   Edad Actividad1 Actividad2 Actividad3 Infarto1 Angina1 Presion1
## 1   52          1          0          0        1       1        1
## 2   66          0          1          0        1       1        1
## 3   56          0          1          0        1       1        1
## 4   57          0          1          0        1       1        1
## 5   42          1          0          0        1       1        1
## 6   62          0          0          0        1       0        1

Modelo general de regresión logísitca

El modelo de regresión logística quedará por lotanto definido de la siguiente manera:

\[ \small \frac{p}{1-p}=\beta_0+\beta_1 \texttt{Edad}+ \beta_2\texttt{Actividad1} +\beta_3\texttt{Actividad2} +\beta_4 \texttt{Actividad3} +\beta_5\texttt{Angina1} +\beta_6\texttt{Presion1} \]

La función \(\tt glm()\)

El model de regresión logística se puede ejecutar con la función \[ \tt glm(modelo,family,data) \]

  • \(\tt modelo\): indica el modelo lineal que queremos aplicar, de manera similar a otros casos que ya hemos visto.
  • \(\tt family\): debemos indicar la familia que utilizaremos en la construcción del modelo. En este caso es \(\tt family=binomial\) (para otras alternativas, las principales son \(\tt gaussian,poisson,Gamma\)).
  • \(\tt data\): Los datos incluidos en data tienen que estar en formato data frame.

Aplicación de la función

mod1 <- glm(Infarto1 ~ Edad+Actividad1+Actividad2+Actividad3+Angina1+Presion1, family = "binomial", valores)


summary(mod1)
## 
## Call:
## glm(formula = Infarto1 ~ Edad + Actividad1 + Actividad2 + Actividad3 + 
##     Angina1 + Presion1, family = "binomial", data = valores)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.21206    2.80797  -0.788 0.430826    
## Edad         0.01706    0.04491   0.380 0.704056    
## Actividad1   2.49278    1.34671   1.851 0.064168 .  
## Actividad2   1.60335    1.19895   1.337 0.181126    
## Actividad3   0.58177    1.29155   0.450 0.652388    
## Angina1     -1.15875    1.28787  -0.900 0.368258    
## Presion1     3.27150    0.95440   3.428 0.000609 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 64.104  on 49  degrees of freedom
## Residual deviance: 41.484  on 43  degrees of freedom
## AIC: 55.484
## 
## Number of Fisher Scoring iterations: 5

Ejectuando el model, se puede ver como la única variable significativa en la Probabilidad de Infarto es la variable \(\tt Presion\).

Con lo cual volveremos a ehecutar el modelo aislado la variable.

mod2 <- glm(Infarto1 ~ Presion1, family = "binomial", valores)


summary(mod2)
## 
## Call:
## glm(formula = Infarto1 ~ Presion1, family = "binomial", data = valores)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7732     0.4935  -1.567 0.117209    
## Presion1      2.6827     0.7284   3.683 0.000231 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 64.104  on 49  degrees of freedom
## Residual deviance: 47.541  on 48  degrees of freedom
## AIC: 51.541
## 
## Number of Fisher Scoring iterations: 4

Utilizando los coeficientes facilitados por la función \(\tt glm()\), el modelo por lo tanto será

\[ \frac{p}{1-p}=\beta_0 + \beta_6·PRESION=-0.7732+2.6827·PRESION \]

Con un \(\texttt{p-value}\) tan bajo que podemos considerar el modelo como ajustado.

\[ p=\frac{\text{exp}(-0.7732+2.6827·PRESION)}{1+\text{exp}(-0.7732+2.6827·PRESION)} \]

ques la probabilidad de \(1\), esto es de éxito (en este caso es infarto Sí).

Probabilidad de infarto con presión baja

La probabilidad de tener infarto con presión baja (PRESION=0) será:

\[ p=\frac{\text{exp}(-0.7732+2.6827·0)}{1+\text{exp}(-0.7732+2.6827·0)}= \frac{\text{exp(-0.7732)}}{1+\text{exp(-0.7732)}}=0.3157873 \]

Probabilidad de infarto con presión alta

La probabilidad de tener infarto con presión alta (PRESION=1) será:

\[ p=\frac{\text{exp}(-0.7732+2.6827)}{1+\text{exp}(-0.7732+2.6827)}=\frac{6.749713}{1+6.749713}=0.870963 \]

Otros Modelos Lineales Generalizados

El modelo \(Logit\) hemos visto que está ligado con las covariables independientes por la función:

\[ g(p)=\text{log}\frac{p}{1-p} \]

El modelo \(Probit\), bastante similar

\[ g(p)=\Phi^{-1}(p) \]

El Modelo de Regresion Poisson

\[ g(p)=\text{log}\mu \]

donde \(\mu\) es la media de las variable dependiente, \(Y\), que sigue la distribución de Poisson.

Modelo \(Probit\) con R

Muy similar a \(Logit\). Solo hay que ejecutar el parámetro correspondiente \(\tt family=binomial(link=probit)\)

Vemos que el ejemplo anterior quedaría así:

mod3 <- glm(Infarto1 ~ Edad+Actividad1+Actividad2+Actividad3+Angina1+Presion1, family = binomial(link=probit), valores)


summary(mod3)
## 
## Call:
## glm(formula = Infarto1 ~ Edad + Actividad1 + Actividad2 + Actividad3 + 
##     Angina1 + Presion1, family = binomial(link = probit), data = valores)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.27810    1.59787  -0.800 0.423781    
## Edad         0.01083    0.02530   0.428 0.668604    
## Actividad1   1.48657    0.75429   1.971 0.048743 *  
## Actividad2   0.88878    0.66098   1.345 0.178742    
## Actividad3   0.35301    0.73562   0.480 0.631312    
## Angina1     -0.75917    0.75229  -1.009 0.312907    
## Presion1     1.93701    0.51233   3.781 0.000156 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 64.104  on 49  degrees of freedom
## Residual deviance: 41.255  on 43  degrees of freedom
## AIC: 55.255
## 
## Number of Fisher Scoring iterations: 6

Comparando \(Logit\) y \(Probit\), por ser la regresión logística de colas un poco más pesadas que la normal, los estimadores que se obtienen con \(Logit\) son algo mayores que los de \(Probit\) (1.6-1.8 veces mayores, aprox.).







  1. \(\text{exp}(x)\) y \(e^x\) son exactamente lo mismo. Ambas expresiones representan la función exponencial natural, donde la constante \(e\) (número de Euler \(\approx 2.718\)) se eleva a la potencia (\(x\)). “Exp” es simplemente la notación abreviada para “exponencial” y se utiliza indistintamente con \(e^x\) en matemáticas, cálculo y programación.↩︎