EL CASO BINARIO
Regresión logística (solo variables categóricas)
Dr. rer. nat. Humberto LLinás Solano
Departamento de Matemáticas y Estadística, Universidad del Norte (Barranquilla, Colombia)14/07/22
Abstract
La teoría mencionada puede revisarse en el volumen 8 de mis notas de clase que aparecen en el siguiente documento: 2.2. Regresión logística y en la referencia: LLinás (2006). En Rpubs:: toc se pueden ver otros documentos de posible interés.
1 Librerías
library(aplore3) #Base de datos para los ejemplos
library(lsm) #Base de datos para ejemplos y estimaciones del Log-verosimilitud
library(tidyverse) #Incluye a dplyr y ggplot22 Introducción
Los métodos de regresión se han convertido en un componente integral de cualquier análisis de datos preocupado por describir la relación entre una variable de respuesta y una o variables más explicativas. Muy a menudo, la variable de resultado es discreta, tomando un valor de dos o más valores posibles. El modelo de regresión logística es el más modelo de regresión de mayor uso frecuente para el análisis de estos datos.
En el documento Rpbus :: Modelos lineales generalizados se explicó que estos modelos hacen parte de los modelos lineales generalizados. Allí se pueden ver los detalles correspondientes. Para conocer con profundidad estos modelos, es importante estudiar los siguientes cuatro tipos de modelos:
Modelo de Bernoulli.
Modelo completo.
Modelo nulo.
Modelo saturado.
En el documento Rpbus :: Modelos completo, nulo y saturado se describieron sus propiedades, con los ejemplos correspondientes y en los documentos Rpbus :: Modelos logísticos (estimaciones), Rpbus :: Modelos logísticos (intervalos de confianza) y Rpbus :: Modelos logísticos (pruebas de hipótesis) se detalló todo lo relacionado con las estimaciones e intervalos de confianza para los parámetros logísticos, así como diferentes pruebas de comparación entre modelos. En este documento, se utilizarán las notaciones utilizados allá, así como los resultados encontrados en los ejemplos aplicados en esos documentos.
3 Datasets
Para las aplicaciones, se utilizaron bases de datos de las librerías aplore3 (creado por Braglia, 2016) y lsm (creado por LLinás, Fábregas y Villalba, 2020).
library(aplore3)
library(lsm)Ambos paquetes incluyen, de manera “no oficial”, todos los conjuntos de datos utilizados en el texto de Hosmer, Lemeshow and Sturdivant (2013). En este link o en este otro se encuentran los nombres de los datasets con los respectivos detalles:
Descripción (description).
Uso (usage).
Formato (format).
Fuente (source) o Referencias (References).
Se resalta el hecho que lsm es un paquete que contiene:
Otros datasets que son de autoría propia.
La función lsm(), la cual nos permite estimar, entre otros, el logaritmo de la función de verosimilitud de los modelos completos, nulo, saturado y logístico. La teoría relacionada con los tres primeros modelos se van a explicar más adelante dentro de este documento. La que se relaciona con lo smodelos logísticos, en otros documentos.
Para los ejemplos, se utilizará chdage, el cual contiene datos recogidos con el fin de estudiar si la edad es un factor influyente en la presencia o no de enfermedades coronarias (CHD). Es un data frame con 100 observaciones. En aplore3 hay cuatro variables y en lsm, tres. A continuación, se describen cada una (entre paréntesis, el nombre de la variable en lsm):
id (ID): código de identificación (1-100).
chd (CHD): Presencia of CHD. En aplore3, los niveles son No, Si. En lsm son 0 (=No) y 1 (=Si).
age (AGE): edad (en años) de los participantes.
agegrp: edad agrupada (1: 20-39, 2: 30-34, 3: 35-39, 4: 40-44, 5: 45-49, 6: 50-54, 7: 55-59, 8: 60-69).
Cargamos el data frame chdage de la librería aplore3:
chdage <- aplore3::chdage
attach(chdage)Las primeras 10 observaciones son:
| ID | age | agegrp | chd |
|---|---|---|---|
| 1 | 20 | 20-39 | No |
| 2 | 23 | 20-39 | No |
| 3 | 24 | 20-39 | No |
| 4 | 25 | 20-39 | No |
| 5 | 25 | 20-39 | Yes |
| 6 | 26 | 20-39 | No |
| 7 | 26 | 20-39 | No |
| 8 | 28 | 20-39 | No |
| 9 | 28 | 20-39 | No |
| 10 | 29 | 20-39 | No |
Observe que, en el datasets aplore3::chdage, la variable chd es binaria, pero no numérica. Al convertirla en integer, se codifica como 1=No y 2=Si. Por esta razón, esta variable se puede codificar como 0=No y 1=Si, de la siguiente manera:
CHD <- as.integer(chd)-14 El caso de variables independientes discretas
Si alguna de las variables independientes son discretas o variables en escala nominal tales como sexo, raza, sexo, grupo de tratamiento, etc., es inapropiado incluirlas en el modelo como si ellas fuesen variables en escala de intervalo. Los números usados para representar los diferentes niveles son identificaciones y no tienen significado numérico. En esta situación, el método de escogencia es usar una colección de variables de diseño (o variables dummy).
Example 4.1 Suponga que una de las variables independientes es raza que ha sido codificado como blanco, negro y otro. En este caso, se necesitan dos variables de diseño, digamos, \(D_1\) y \(D_2\) (véase la Tabla 4.1).
| 1 | 2 | 3 |
|---|---|---|
| Raza | \(D_1\) | \(D_2\) |
| Blanca | \(0\) | \(0\) |
| Negra | \(1\) | \(0\) |
| Otra | \(0\) | \(1\) |
Remark. Tenemos:
En general, si una variable escala nominal tiene \(k\) posibles valores, entonces se necesitan \(k-1\) variables de diseño.
Para ilustrar la notación usada para las variables de diseño en estas notas, suponga que la \(j\)-ésima variable independiente \(x_j\) tiene \(k_j\) niveles. Las \(k_j-1\) variables de diseño serán denotadas por \(D_{jl}\) y los coeficientes para estas variables de diseño como \(\beta_{jl}\). Por consiguiente, el logit para un modelo probabilístico con \(k\) variables y la \(j\)-ésima variable discreta será:
\[\begin{eqnarray} g(x):= \mbox{Logit}(p_j)= \delta \;+\; \beta_1 \,x_{1} \;+\;\cdots \;+\; \sum\limits_{m=1}^{k_j-1}\beta_{jm} D_{jm} \;+\; \beta_K \,x_{K} \tag{4.1} \end{eqnarray}\]
5 Variable independiente dicotómica
5.0.1 Descripción del caso
Supongamos que tenemos, por ejemplo, un estudio de casos y control, donde codificamos con:
\(y=1\), si la persona pertenece al grupo de los casos (es decir, si tiene la enfermedad).
\(y=0\), si la persona pertenece al grupo control (es decir, si no tiene la enfermedad).
Además, consideraremos el caso donde el modelo logístico contiene sólo una variable independiente y que ésta es nominal y dicotómica (es decir, medida en dos niveles). Asumamos que la variable independiente \(x\) está codificada como:
\(x=0\), si la persona no está expuesta a un factor de riesgo.
\(x=1\), si la persona está expuesta a un factor de riesgo.
Para simplificar la notación, haremos uso de:
\[\pi(x):= P(Y=1/X=x)\]
para representar la probabilidad condicional de \(Y=1\) dado \(X=x\) cuando se utiliza la regresión logística. En este caso, por el supuesto 4 (ver Hipótesis 5.2 en la sección 5 del documento Rpbus :: Modelos logísticos (estimaciones)), las dos ecuaciones del modelo serán:
\[\begin{eqnarray*} g(0) &:=& \ln\left(\frac{\pi(0)}{1-\pi(0)}\right) = \delta \;+\; \beta \cdot 0 \;=\; \delta\\ &&\\ g(1) &:=& \ln\left(\frac{\pi(1)}{1-\pi(1)}\right) = \delta \;+\; \beta \cdot 1 \;=\; \delta \;+\; \beta \end{eqnarray*}\]
Como ya se explicó en el Teorema 6.1 (sección 6) del documento Rpbus :: Modelos logísticos (estimaciones), viene dada por
\[\pi(x) \;= \; \frac{e^{g(x)}}{1+ e^{g(x)}}, \qquad x\in\{0,1\}\] Los posibles valores de las probabilidades logísticas pueden ser convenientemente organizarse en una tabla de \(2\times 2\) como, se muestra en la Tabla @ref(tab:Tabla-3.1-Hosmer-pag-49}. En ella vemos las formulas de \(\pi(x)\) para cada valor de \(x\) y de \(y\):
| 1 | 2 | 3 |
|---|---|---|
| Valores de \(x\) y \(y\) | \(x=1\) (expuesta) | \(x=0\) (no expuesta) |
| \(y=1\) (caso) | \(\pi(1) \;= \; \frac{e^{\delta + \beta}}{1+ e^{\delta + \beta}}\) | \(\pi(0) \;= \; \frac{e^{\delta}}{1+ e^{\delta}}\) |
| \(y=0\) (control) | \(1-\pi(1) \;= \; \frac{1}{1+ e^{\delta + \beta}}\) | \(1-\pi(0) \;= \; \frac{1}{1+ e^{\delta}}\) |
| Total | \(1\) | \(1\) |
Remark. Debemos enfatizar que el primer paso, para interpretar el efecto de la variable independiente en un modelo, es expresar la diferencia de logit deseada en términos del modelo. Observe que, en este caso, esta diferencia es igual a \(\beta\):
\[g(1) - g(0) \;= \; (\delta + \beta_1) - \delta \;= \; \beta\]
Para interpretar este resultado necesitamos introducir y discutir la razón odds (véase la sección siguiente), conceptos que fueron descritos con detalles en el documento Rpbus :: Modelos logísticos (estimaciones).
5.0.2 Odds y razón odds
Definition 5.1 (Odds) Los odds de los resultados cuando están presente individuos con \(x=1\) está definido como
\[O(1)=\frac{\pi(1)}{1-\pi(1)}\]
Similarmente, los odds de los resultados cuando están presente individuos con \(x=0\) está definido como
\[O(0)=\frac{\pi(0)}{1-\pi(0)}\]
Definition 5.2 (Razón odds) La se define como el cociente entre el odds para \(x=1\) y el odds para \(x=0\) y está dada por la ecuación:
\[\begin{eqnarray} OR \;= \; \frac{O(1)}{O(0)} \;= \; \frac{\pi(1)/[1-\pi(1)]}{\pi(0)/[1-\pi(0)]} \end{eqnarray}\]
Theorem 5.1 (Razón odds con exponencial de la pendiente) Teniendo en cuenta las expresiones la Tabla 5.1 y la Definición 5.2, obtenemos:
\[ OR \;=\; \frac{\left(\frac{e^{\delta + \beta}}{1+ e^{\delta + \beta}}\right) \Big/ \left(\frac{1}{1+ e^{\delta + \beta}}\right)}{\left( \frac{e^{\delta }}{1+ e^{\delta}}\right) \Big/ \left( \frac{1}{1+ e^{\delta }}\right)}\;=\; \frac{e^{\delta + \beta}}{e^{\delta}} \;=\; e^{(\delta + \beta) - \delta} \;=\; e^{\beta} \]
Remark. Tenemos:
Para la regresión logística con una variable independiente codificada con 0 y 1, la relación entre la razón odds y el coeficiente de regresión es como se muestra en el Teorema 5.1.
Si \(OR=k\), una interpretación sería cualquiera de las dos propuestas de abajo:
La razón entre la presencia (casos, \(y=1\)) versus la no presencia de la enfermedad (controles, \(y=0\)) es \(k\) veces mayor en las personas expuestas a un factor de riesgo (población \(x=1\)) en comparación a las personas no expuestas (población \(x=0\)).
La presencia de la enfermedad (casos, \(y=1\)) es \(k\) veces más probable que ocurra en las personas expuestas al factor (población \(x=1\)) más que las personas no expuestas (población \(x=0\)).
Example 5.1 Suponga que \(Y\) representa la presencia (\(y=1\)) o ausencia de cáncer de hígado (\(y=0\))y \(X\) denota si una persona es fumadora (\(x=1\)) o no (\(x=0\)). Entonces un valor \(\widehat{OR}=2\) estima que el cáncer de hígado es dos veces más probable que ocurra entre los fumadores que entre los no fumadores en la población de estudio.
Example 5.2 Suponga que \(Y\) representa la presencia (\(y=1\)) o ausencia de infarto (\(y=0\)) y \(X\) denota si una persona está comprometida a realizar ejercicios físicos de manera extenuante (\(x=1\)) o no (\(x=0\)). Si el odds estimado es \(\widehat{OR}=0.5\), entonces la ocurrencia de infarto es 0.5 más probable que ocurra entre los que hacen ejercicio que entre aquéllos que no lo hacen en la población de estudio. O se puede decir, que la ocurrencia de infarto es 2 veces más probable que ocurra enrte los que no hacen ejercicio que entre los que lo hacen en la población de estudio.
6 Ejemplo 1: Enunciado
Considere los datos del archivo chdage.
Crear una nueva variable, AGED, que toma valores 1 si la edad del sujeto es mayor o igual que 55 y 0, de otro modo.
Construir una tabla de contingencia (es decir, un tabulación cruzada) entre la variable edad dicotomizada AGED con la variable de respuesta CHD.
Halle el logaritmo de la función de verosimilitud \({\cal L}(p)\) en el modelo saturado.
Halle la función de verosimilitud \(L(p)\) en el modelo saturado.
Escriba las ecuaciones del modelo logístico estimado (solo indicando los estimadores).
Estime los parámetros logísticos.
Estime los respectivos errores estándares de los estimadores de los parámetros logísticos con un sofware estadístico.
Estime los respectivos errores estándares de los estimadores de los parámetros logísticos de manera manual.
Con los resultado encontrados en (f), escriba las ecuaciones estimadas de regresión logística.
Encuentre \(\hat{\pi}(0)\) y \(\hat{\pi}(1)\).
Con los resultados encontrados en (j), calcule \({\cal L}(\tilde{p})\).
Con los resultados encontrados en (j), calcule \(L(\tilde{p})\).
Encuentre la estimación de la razón odds.
Construya un intervalo del \(95\%\) para la razón odds.
7 Ejemplo 1: Solución
7.0.1 Solución parte (a)
Vamos a crear una nueva variable, AGED, que toma valores 1 si la edad del sujeto es mayor o igual que 55 y 0 de otro modo.
datos <- lsm::chdage
attach(datos)
datos$AGED <- ifelse(AGE<55, 0, 1)
AGED <- datos$AGED7.0.2 Solución parte (b)
El resultado de cruzar la variable edad dicotomizada AGED con la variable de respuesta CHD se presenta en la Tabla 7.1:
table(CHD, AGED)| 1 | 2 | 3 | 4 |
|---|---|---|---|
| CHD (\(y\)) | AGED \(\geq 55\) (\(x=1\)) | AGED < \(55\) (\(x=0\)) | Total |
| Presente (\(y=1\)) | \(21\) | \(22\) | \(43\) |
| Ausente (\(y=0\)) | \(6\) | \(51\) | \(57\) |
| Total | \(27\) | \(73\) | \(100\) |
7.0.3 Solución parte (c)
El modelo saturado tiene \(J=2\) poblaciones, que son los grupos \(x=0\) y \(x=1\). Entonces, por el teorema 4.1 que aparece en la sección 4 del documento Rpbus :: Modelos logísticos (estimaciones), se tiene que:
\[\begin{eqnarray*} {\cal L}(p) &=& \sum^2_{j=1}[z_j\ln p_j + (n_j- z_j)\ln (1-p_j)] \\ &=& \big\{z_1\ln \pi(0) \;+\; (n_1- z_1)\ln \big(1-\pi(0)\big)\big\} \;+\; \big\{z_2\ln \pi(1) \;+\; (n_2- z_2)\ln \big(1-\pi(1)\big)\big\}\\ &&\\ &=& 22\ln \pi(0) \;+\; 51\ln \big(1-\pi(0)\big) \;+\; 21\ln \pi(1) \;+\; 6\ln \big(1-\pi(1)\big) \end{eqnarray*}\]
7.0.4 Solución parte (d)
Sabiendo que
\[{\cal L}(p) \;=\; \ln \big(L(p)\big) \]
entonces, por el inciso (c), la función de verosimilitud en el modelo saturado es:
\[L(p) \;= \; \pi(1)^{21} \cdot [1-\pi(1)]^6 \cdot \pi(0)^{22} \cdot [1-\pi(0)]^{51}\]
7.0.5 Solución parte (e)
Las ecuaciones estimadas son:
\[\begin{eqnarray*} \hat{g}(0) &:=& \ln\left(\frac{\hat{\pi}(0)}{1 - \hat{\pi}(0)}\right) = \hat{\delta} \;+\; \hat{\beta}_{AGED=0} \cdot 0 \;=\; \hat{\delta}\\ &&\\ \hat{g}(1) &:=& \ln\left(\frac{\hat{\pi}(1)}{1 - \hat{\pi}(1)}\right) = \hat{\delta} \;+\; \hat{\beta}_{AGED=1} \cdot 1 \;=\; \hat{\delta} \;+\; \hat{\beta}_{AGED=1} \end{eqnarray*}\]
7.0.6 Solución parte (f)
Con el paquete lsm, podemos obtener las estimaciones solicitadas (ver código y recuadro rojo en la figura 7.1):
\[\hat{\delta}= -0.8408, \qquad \hat{\beta}_{AGED=1}= 2.0935 \]
datos <- lsm::chdage
attach(datos)
datos$AGED <- ifelse(AGE<55, 0, 1)
AGED <- datos$AGED
modelo <- lsm(CHD ~ AGED, data=datos)
modelo
modelo$coefficients # Vector de parámetros
delta <- modelo$coefficients[1] # Intercepto
delta
beta <- modelo$coefficients[2] # Pendiente
beta
Figure 7.1: Estimaciones de los parámetros logísticos. Fuente: Elaboración propia.
Se resalta el hecho que debe cumplirse:
\[\hat{\beta}_{AGED=0} \;+\; \hat{\beta}_{AGED=1}\;=\; 0\]
Por esta razón,
\[\hat{\beta}_{AGED=0} \;=\; - \hat{\beta}_{AGED=1}\;=\; - 2.0935 \]
7.0.7 Solución parte (g)
Con el paquete lsm, podemos obtener las estimaciones solicitadas (ver código y recuadro rojo en la figura 7.2):
\[\hat{S}_{\hat{\delta}} = 0.2551, \qquad \hat{S}_{\hat{\beta}_{AGED=1}} = 0.5285 \]
datos <- lsm::chdage
attach(datos)
datos$AGED <- ifelse(AGE<55, 0, 1)
AGED <- datos$AGED
modelo <- lsm(CHD ~ AGED, data=datos)
modelo
modelo$Std.Error
Figure 7.2: Estimaciones de los errores estándares en el modelo logístico. Fuente: Elaboración propia.
7.0.8 Solución parte (h)
Con base en las frecuencias de las celdas de la Tabla 7.1, se resalta el hecho de que:
- La estimación del error estándar de \(\hat{S}_{\hat{\delta}}\) se puede calcular con la raiz cuadrada de la suma de los inversos multiplicativos de las frecuencias asociadas al grupo \(x=0\):
\[\hat{S}_{\hat{\delta}} \;=\; \sqrt{\frac{1}{22} \;+\; \frac{1}{51}} \;=\; 0.2551\] 2. La estimación del error estándar de \(\hat{S}_{\hat{\beta}_{AGED=1}}\) se puede calcular con la raiz cuadrada de la suma de los inversos multiplicativos de las frecuencias asociadas tanto al grupo \(x=0\) como al grupo \(x=1\):
\[\hat{S}_{\hat{\beta}_{AGED=1}} \;=\; \sqrt{\frac{1}{21} \;+\; \frac{1}{6} \;+\; \frac{1}{22} \;+\; \frac{1}{51} } \;=\; 0.5285\]
SE_delta <- sqrt(1/22 + 1/51)
SE_delta
SE_beta <- sqrt(1/21 + 1/6 + 1/22 + 1/51)
SE_beta7.0.9 Solución parte (i)
Teniendo en cuenta los resultados encontrados en (e) y (f), las ecuaciones estimadas de regresión logística son:
\[\begin{eqnarray*} \hat{g}(0) &:=& \ln\left(\frac{\hat{\pi}(0)}{1 - \hat{\pi}(0)}\right) \;=\; \hat{\delta} \;=\; -0.8408\\ &&\\ \hat{g}(1) &:=& \ln\left(\frac{\hat{\pi}(1)}{1 - \hat{\pi}(1)}\right) \;=\; \hat{\delta} \;+\; \hat{\beta}_{AGED=1} \;=\; -0.8408 \;+\; 2.0935 \;=\; 1.2528 \end{eqnarray*}\]
datos <- lsm::chdage
attach(datos)
datos$AGED <- ifelse(AGE<55, 0, 1)
AGED <- datos$AGED
modelo <- lsm(CHD ~ AGED, data=datos)
delta <- modelo$coefficients[1] # Intercepto
beta <- modelo$coefficients[2] # Pendiente
g_0 <- delta
g_0
g_1 <- delta + beta
g_17.0.10 Solución parte (j)
Por lo explicado en la sección 5.0.1:
\[\begin{eqnarray*} \hat{\pi}(0) &=& \frac{e^{\hat{g}(0)}}{1+ e^{\hat{g}(0)}} \;=\; \frac{e^{-0.8408}}{1+ e^{-0.8408}} \;=\; 0.3014\\ &&\\ \hat{\pi}(1) &=& \frac{e^{\hat{g}(1)}}{1+ e^{\hat{g}(1)}} \;=\; \frac{e^{1.2528}}{1+ e^{1.2528}}\;=\; 0.7778\\ \end{eqnarray*}\]
Con el paquete lsm:
pi_0 <- exp(g_0)/(1+ exp(g_0))
pi_0
pi_1 <- exp(g_1)/(1+ exp(g_1))
pi_1
unique(modelo$p_hat)7.0.11 Solución parte (k)
\[ {\cal L}(\tilde{p}) \;=\; 22\ln \hat{\pi}(0) \;+\; 51\ln \big(1-\hat{\pi}(0)\big) \;+\; 21\ln \hat{\pi}(1) \;+\; 6\ln \big(1-\hat{\pi}(1)\big) \;=\; -58.9796\]
#1. Primera forma (manual):
Log_Lik_Sat <- 22*log(pi_0) + 51*log(1-pi_0) + 21*log(pi_1) + 6*log(1-pi_1)
Log_Lik_Sat
#2. Segunda forma (con lsm):
datos <- lsm::chdage
attach(datos)
datos$AGED <- ifelse(AGE<55, 0, 1)
AGED <- datos$AGED
modelo <- lsm(CHD ~ AGED, data=datos)
Log_Lik_Sat <- modelo$Log_Lik_Saturate
Log_Lik_Sat7.0.12 Solución parte (l)
Sabemos que
\[{\cal L}(\tilde{p}) \;=\; \ln \big(L(\tilde{p})\big) \]
Entonces, la función de verosimilitud en el modelo saturado es:
\[L(\tilde{p}) \;=\; e^{{\cal L}(\tilde{p})} \;=\; 0\]
Otra forma de calcular \(L(\tilde{p})\) es utilizando el inciso (c):
\[L(\tilde{p}) \;= \; \hat{\pi}(0)^{22} \cdot [1-\hat{\pi}(0)]^{51} \cdot \hat{\pi}(1)^{21} \cdot [1-\hat{\pi}(1)]^6 \;=\; 0\]
#1. Primera forma (manual):
Lik_Sat <- exp(Log_Lik_Sat)
Lik_Sat
#2. Segunda forma (manual)
Lik_Sat <- (pi_0)^(22)*(1-pi_0)^(51)*(pi_1)^(21)*(1-pi_1)^6
Lik_Sat 7.0.13 Solución parte (m)
Según el Teorema 5.1, la estimación de la razón de odds es
\[\widehat{OR}\; =\; e^{2.0935} \; =\; 8.1136\]
Este valor puede obtenerse directamente de la Tabla 7.1 o con el paquete lsm, de la siguiente manera (ver código y recuadro rojo en la figura 7.3):
\[\widehat{OR} \;= \; \frac{21/6}{22/51} \;= \; 8.1136\]
datos <- lsm::chdage
attach(datos)
datos$AGED <- ifelse(AGE<55, 0, 1)
AGED <- datos$AGED
modelo <- lsm(CHD ~ AGED, data=datos)
modelo$ExpB[2]
Figure 7.3: Estimación del OR en el modelo logístico. Fuente: Elaboración propia.
7.0.14 Solución parte (n)
Por el teorema 14.7, de la sección 14 del documento Rpbus :: Modelos logísticos (intervalos de confianza), un intervalo de confianza del \((1-\alpha)100\%\) para \(OR\) es:
\[\exp\left\{\hat{\beta}_{AGED=1} \; - \; Z_{\alpha/2} \, \hat{S}_{\hat{\beta}_{AGED=1}} \right\}\quad < \quad OR \quad < \exp\left\{\hat{\beta}_{AGED=1} \; + \; Z_{\alpha/2} \, \hat{S}_{\hat{\beta}_{AGED=1}} \right\}\]
Recuerde que \(\hat{S}_{\hat{\beta}_{AGED=1}}\) es el error estándar del estimador \(\hat{\beta}_{AGED=1}\). Teniendo en cuenta lo anterior, un intervalo de confianza del \(95\%\) para la razón odds poblacional es:
\[\mbox{2.87956} \;<\; \mbox{OR} \;<\; \mbox{22.86148}\]
modelo <- lsm(CHD~AGED, data=datos)
beta <- modelo$coefficients[2]
ES_beta <- modelo$Std.Error[2]
alfa <- 0.05
Z <- qnorm(1-alfa/2)
CI_lower_beta <- beta - Z*ES_beta
CI_upper_beta <- beta + Z*ES_beta
CI_beta <- cbind(CI_lower_beta, CI_upper_beta)
#CI_beta
CI_lower <- exp(CI_beta[1])
CI_upper <- exp(CI_beta[2])
CI_OR <- cbind(CI_lower, CI_upper)
CI_OREste intervalo es sesgado a la derecha porque excede el 1. Entonces con un grado de confianza del \(95\%\), el CHD ocurre entre aquellas personas con edad mayor o igual que 55 años en la población de estudio por lo menos 2.9 veces o a lo más 22.9 veces más probable que aquéllos por debajo de 55.
8 Variable independiente policotómica
En esta sección supondremos que, en vez de dos categorías, la variable independiente tiene \(k>2\) valores diferentes. Analizaremos la situación a través del siguiente ejemplo.
8.0.1 Ejemplo: Enunciado
Suponga un caso hipotético. En un estudio de enfermedades coronarias se tienen en cuenta solo a las variables CHD (\(Y\)) (codificada como \(y=1\), si tiene la enfermedad y \(y=0\), en caso contrario) y Raza (la cual está codificada en cuatro niveles). La tabulación cruzada de Raza por CHD produce los datos que aparecen en la Tabla 8.1.
CHD <- c(rep(1, 5), rep(0, 20), rep(1, 15), rep(0, 10), rep(1, 20), rep(0, 10), rep(1, 10), rep(0, 10))
Raza <- c(rep("Blanca", 25), rep("Hispana", 25), rep("Negra", 30), rep("Otra", 20))
chdraza <- data.frame(CHD,Raza)
attach(chdraza)
head(chdraza)
table(CHD, Raza)| 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|
| No. | CHD (\(y\)) | Blanca | Hispana | Negra | Otra | Total |
| \(1.\) | Presente (\(y=1\)) | \(5\) | \(15\) | \(20\) | \(10\) | \(50\) |
| \(2.\) | Ausente (\(y=0\)) | \(20\) | \(10\) | \(10\) | \(10\) | \(50\) |
| \(3.\) | Total | \(25\) | \(25\) | \(30\) | \(20\) | \(100\) |
| ESTIMACIONES: | ||||||
| 4. Las relacionadas con el intercepto son \(\hat{\delta}=-1.3863\) y \(\hat{S}_{\hat{\delta}}=0.500\) | ||||||
| \(5.\) | \(\hat{\beta}\) | \(-5.2575\) | \(1.7918\) | \(2.0794\) | \(1.3863\) | |
| \(6.\) | \(\hat{S}_{\hat{\beta}}\) | \(0.6455\) | \(0.6325\) | \(0.6708\) | ||
| \(7.\) | \(\widehat{OR}\) | \(1\) | \(6\) | \(8\) | \(4\) | |
| \(8.\) | \(\ln(\widehat{OR}) = \hat{\beta}\) | \(0\) | \(1.7918\) | \(2.0794\) | \(1.3863\) | |
| \(9.\) | IC \(95\%\) para OR | \((1.7; 21.3)\) | \((2.3; 27.6)\) | \((1.1; 14.9)\) | ||
| Fuente. Elaboración propia. | ||||||
| Notaciones. 1 \(\widehat{OR}:\) Razón odds estimada; 2 \(OR:\) Razón odds; 3 IC \(95\%\): Intervalo del \(95\%\) de confianza. | ||||||
| Intercepto. a \(\hat{\delta}\): su estimación; b \(\widehat{SE}_{\hat{\delta}}\): estimación de su error estándar. | ||||||
| Pendiente. * \(\hat{\beta}\): su estimación; † \(\widehat{SE}_{\hat{\beta}}\): estimación de su error estándar. | ||||||
En esa misma tabla se dan los totales (fila 3 y columna 7), las estimaciones de los parámetros logísticos y los errores de sus respectivos estimadores (filas 4, 5 y 6), la razón odds estimada para cada raza usando Blanco como grupo de referencia (fila 7), el intervalo de confianza para cada razón odds OR (fila 8) y el logaritmo de cada OR estimado (fila 9). A continuación, de manera muy general se darán las justificaciones correspondientes.
8.0.2 Estimaciones de \(\delta\) y \(\beta\)
las estimaciones de los parámetros logísticos se pueden obtener, por ejemplo, con glm: (ver recuadro rojo en la figura 8.1):
\[\hat{\delta}= -1.3863, \qquad \widehat{\beta}_{hispana}= 1.7918, \qquad \widehat{\beta}_{negra}= 2.0794, \qquad \widehat{\beta}_{otra}= 1.3863\]
modelo <- glm(CHD~Raza, family=binomial(link = "logit"), data=chdraza)
summary(modelo)
coefficients(modelo)
Figure 8.1: Estimación de parámetros logísticos. Fuente: Elaboración propia.
Se resalta el hecho que debe cumplirse:
\[\widehat{\beta}_{blanca}\;+\; \widehat{\beta}_{hispana} \;+\; \widehat{\beta}_{negra} \;+\; \widehat{\beta}_{otra} \;=\; 0\]
Por esta razón,
\[\widehat{\beta}_{blanca} \;=\; - \left(\widehat{\beta}_{hispana} \;+\; \widehat{\beta}_{negra} \;+\; \widehat{\beta}_{otra}\right)\;=\; - 5.2575 \]
8.0.3 Estimaciones de los errores estándares con lsm
las estimaciones de los errores estándares de los estimadores de los parámetros logísticos se pueden obtener, por ejemplo, con glm: (ver recuadro rojo en la figura 8.2).
\[\hat{S}_{\hat{\delta}} = 0.5, \qquad \hat{S}_{\widehat{\beta}_{hispana}} = 0.6455, \qquad \hat{S}_{\widehat{\beta}_{negra}} = 0.6325, \qquad \hat{S}_{\widehat{\beta}_{otra}} = 0.6708 \] Observe que los errores estándares son las raices cuadradas de los elementos diagonales de la matriz de varianzas y covarianzas.
modelo <- glm(CHD~Raza, family=binomial(link = "logit"), data=chdraza)
summary(modelo)
VarCov <- vcov(modelo) #Matriz de varianzas y covarianzas
sqrt(VarCov) #Raiz de esta matriz
Figure 8.2: Estimación de los errores estándares. Fuente: Elaboración propia.
8.0.4 Estimaciones de los errores estándares manual
Con base en las frecuencias de las celdas de la Tabla 8.1, se resalta el hecho de que:
- La estimación del error estándar de \(\hat{S}_{\hat{\delta}}\) se puede calcular con la raiz cuadrada de la suma de los inversos multiplicativos de las frecuencias asociadas al grupo de referencia (blanco):
\[S_{\hat{\delta}}\;= \;\sqrt{\frac{1}{5} + \frac{1}{20}} \;= \; 0.5\]
SE_delta <- sqrt(1/5 + 1/20)
SE_delta- La estimación del error estándar de \(\hat{S}_{\hat{\beta}_{hispana}}\) se puede calcular con la raiz cuadrada de la suma de los inversos multiplicativos de las frecuencias asociadas al grupo de referencia (blanco) y a hispano:
\[\hat{S}_{\hat{\beta}_{hispana}}\;= \;\sqrt{\frac{1}{5} + \frac{1}{20} +\frac{1}{15} +\frac{1}{10}} \;= \; 0.6455\]
SE_beta_hispano <- sqrt(1/5 + 1/20 + 1/15 + 1/10)
SE_beta_hispano- La estimación del error estándar de \(\hat{S}_{\hat{\beta}_{negra}}\) se puede calcular con la raiz cuadrada de la suma de los inversos multiplicativos de las frecuencias asociadas al grupo de referencia (blanco) y a negro:
\[\hat{S}_{\hat{\beta}_{negra}}\;= \;\sqrt{\frac{1}{5} + \frac{1}{20} +\frac{1}{20} +\frac{1}{10}} \;= \; 0.6325\]
SE_beta_negro <- sqrt(1/5 + 1/20 + 1/15 + 1/10)
SE_beta_negro- La estimación del error estándar de \(\hat{S}_{\hat{\beta}_{otra}}\) se puede calcular con la raiz cuadrada de la suma de los inversos multiplicativos de las frecuencias asociadas al grupo de referencia (blanco) y a otro:
\[\hat{S}_{\hat{\beta}_{otra}}\;= \;\sqrt{\frac{1}{5} + \frac{1}{20} +\frac{1}{10} +\frac{1}{10}} \;= \; 0.6708\]
SE_beta_otro <- sqrt(1/5 + 1/20 + 1/10 + 1/10)
SE_beta_otro8.0.5 Estimaciones de los OR
- Para los de raza hispana, la razón odds estimada es
\[\widehat{OR}(\mbox{hispana vs blanca}) \;=\; \frac{(15)(20)}{(10)(5)} \;= \; 6\] 2. Para los de raza negra, la razón odds estimada es
\[\widehat{OR}(\mbox{negra vs blanca}) \;=\; \frac{(20)(20)}{(10)(5)} \;= \; 8\]
- Para los de otra raza, la razón odds estimada es
\[\widehat{OR}(\mbox{otra vs blanca}) \;=\; \frac{(10)(20)}{(10)(5)} \;= \; 4\]
La referencia de grupo es indicada por un valor de 1 para la razón odds. Estas mismas estimaciones para las razones odds, se pueden obtener, por ejemplo, con glm:
modelo <- glm(CHD~Raza, family=binomial(link = "logit"), data=chdraza)
exp(coefficients(modelo))8.0.6 Intervalos de confianza
Por el teorema 14.2, de la sección 14 del documento Rpbus :: Modelos logísticos (intervalos de confianza), un intervalo de confianza del \((1-\alpha)100\%\) para los coeficientes son de la forma:
\[\widehat{\beta}_j \;\pm \; Z_{\alpha/2} S_{\widehat{\beta}_j}\] Por el teorema 14.7, de la sección 14 de ese mismo documento, un intervalo de confianza del \((1-\alpha)100\%\) para \(OR\) es:
\[\exp\{\widehat{\beta}_j \;\pm \; Z_{\alpha/2} S_{\widehat{\beta}_j}\}\]
Teniendo en cuenta lo anterior, un intervalo de confianza del 95% para la pendiente \(\beta_{hispana}\) y para la razón odds \(OR(\mbox{hispana vs blanca})\) son:
\[ 0.5266\;<\; \beta_{hispana} \; <\; 3.0569,\qquad 1.6932 \;<\; OR(\mbox{hispana vs blanca}) \; <\; 21.2618\]
modelo <- glm(CHD~Raza, family=binomial(link = "logit"), data=chdraza)
beta <- modelo$coefficients[2]
beta
ES <- sqrt(vcov(modelo))[2,2]
ES
alfa <- 0.05
Z <- qnorm(1-alfa/2)
CI_lower_beta <- beta - Z*ES
CI_upper_beta <- beta + Z*ES
CI_beta <- cbind(CI_lower_beta, CI_upper_beta)
CI_beta
CI_lower <- exp(CI_beta[1])
CI_upper <- exp(CI_beta[2])
CI_OR <- cbind(CI_lower, CI_upper)
CI_ORSe deja al lector el cálculo de los otros intervalos de confianza presentados en la Tabla 8.1.
8.0.7 Pruebas de comparación
Se deja al lector las pruebas de comparación del modelo logístico con los modelos nulo, saturado y completo. Para ello, se puede apoyar en las salidas que se obtienen con glm::summary() (ver recuadro rojo en la figura 8.3):
modelo <- glm(CHD~Raza, family=binomial(link = "logit"), data=chdraza)
summary(modelo)
Figure 8.3: Salida con glm::summary(). Fuente: Elaboración propia.
9 Ejercicios
Para la solución de los siguientes ejercicios, téngase en cuenta los siguientes comentarios:
Todos los datos mencionados aparecen en los paquetes mencionados en este documento.
Siempre debe detallar el análisis del conjunto de datos (con las variables especificadas) basado en lo explicado en este documento.
Verifique cómo se obtienen las estimaciones correspondientes, los logaritmos de las funciones de máxima verosimilitud, estimaciones e intervalos de confianza para \(p_j\), ODDS, razones ODDS, intercepto, pendiente, etc.
9.0.1 Ejercicios 1 a 3
Demuestre siguientes teoremas del documento Rpbus :: Modelos logísticos (pruebas de hipótesis): (a) 5.1; (b) 6.2; (c) 7.1; (d) 7.2; (e) 9.1.
Haga un listado de los paquetes de R que estimen el logaritmo de la función de máxima verosimilitud en los modelos completo y nulo.
Haga un listado de los paquetes de R que, en el caso binario, estimen el logaritmo de la función de máxima verosimilitud en los modelos saturado y logístico.
9.0.2 Ejercicios 4 a 5
- Los datos ICU corresponden a una muestra de 200 sujetos que hicieron parte de un estudio de supervivencia de pacientes que fueron remitidos a una unidad de cuidados intensivos (intensive care unit - ICU). La meta principal de este estudio fue desarrollar un modelo de regresión logística para predecir la probabilidad de supervivencia de estos pacientes en el hospital y estudiar los factores de riesgos asociados con el índice de mortalidad ICU. Crear una nueva variable, AGED, que toma el valor 1 si la edad es mayor o igual que 60 y el valor 0, de otro modo. Repita todos los análisis realizados en este documento, pero considerando AGED como variable independiente y STA como variable dependiente.
- Considere los datos ICU. Repita los inciso (a), (b), (d), (e) y (g) del ejercicio 4 del documento Rpbus :: Modelos logísticos (pruebas de hipótesis) utilizando la variable TYP como variable independiente (en vez de AGE) y STA como variable dependiente.
9.0.3 Ejercicio 6 a 7
- Considere los datos ICU. Use STA como variable dependiente y CPR como variable independiente.
Demuestre que el valor del logaritmo de las razones odds obtenido de la tabla de contingencia de STA contra CPR es idéntico a la pendiente estimada del modelo de regresión logístico de STA sobre CPR.
Verifique que el error estándar estimado de la pendiente estimada para CPR obtenido en Statgraphics es idéntico a la raiz cuadrada de la suma de los inversos de las frecuencias de cada celda de la tabla de contingencia de STA contra CPR.
¿Qué aspecto relacionado con la codificación de la variable CPR hace que los cálculos para los dos métodos en (a) y (b) sean equivalente.
Obtenga un intervalo del \(95\%\) de confianza para la razón odds.
Para propósitos de ilustración, use una transformación de datos para recodificar (sólo en este inciso) la variable CPR como sigue: 4= no y 2=si. Desarrolle la regresión logística de STA sobre CPR (recodificada). Demuestre cómo los cálculos de la diferencia logit de CPR=si versus CPT=no es equivalente al valor del logaritmo de la razón odds obtenido en (a). Use los resultados de la regresión logística con el fin de obtener un intervalo del \(95\%\) de confianza para la razón odds y verificar que ellos son los mismos límites obtenidos en (d).
- Considere los datos ICU. Use STA como variable dependiente y RACE como variable independiente.
La variable RACE está codificada en tres niveles. Prepare una tabla mostrando la codificación de dos variables de diseño usando como grupo de referencia el valor RACE=1 (White).
Muestre que los logaritmos de las razones odds obtenido de la tabla de contingencia de STA contra RACE, usando RACE=1 como grupo de referencia, son idénticos a las pendientes estimadas para las dos variables de diseño del modelo de regresión logístico de STA sobre RACE.
Verifique que los errores estándar estimados de la pendiente estimada para las dos variables de diseño del modelo son idénticos a la raiz cuadrada de la suma de los inversos de las frecuencias de cada celda de la tabla de contingencia de STA contra RACE que se utilizó para calcular las razones odds.
Obtenga un intervalo del \(95\%\) de confianza para las razones odds.
9.0.4 Ejercicios 8 a 9
- Considere los datos ICU. Use STA como variable dependiente y nuevamente RACE como covariable.
La variable RACE está codificada en tres niveles. Crear variables de diseño para RACE usando como grupo de referencia el valor RACE=1 (White).
Desarrolle la regresión logística de STA contra RACE.
Muestre que las diferencias logit estimadas de RACE=2 versus RACE=1 y RACE=3 versus RACE=1 son equivalentes a los valores del logaritmo de la razón odds obtenidos en el problema 7b.
Use los resultados de la regresión logística para encontrar el \(95\%\) de confianza para las razones odds y verifique que son los mismos límites a los encontrados en el problema 7d.
Halle la matriz de covarianza estimada y obtenga las varianzas de las diferencias logit estimadas de RACE=2 versus RACE=1 y RACE=3 versus RACE=1. Encuentre los errores estándar correspondientes.
Use los resultados de la regresión logística para obtener un intervalo del \(95\%\) de confianza para las razones odds y verificar que ellos son los mismos límites obtenidos en 7d.
- Considere la variable AGE en los datos ICU. Use STA como variable dependiente.
Prepare una tabla mostrando la codificación de tres variables de diseño basadas en los cuartiles empíricos de AGE usando el primer cuartil como referencia de grupo.
Ajuste la regresión logística de STA sobre AGE (recodificada).
Grafique las tres pendientes estimadas versus el punto medio de la edad cuartil respectiva y grafique como un cuarto punto un valor de cero en el punto medio del primer cuartil de la edad. ¿Sugiere esta gráfica que el logit es lineal en la edad?
9.0.5 Ejercicio 10
- Use los datos ICU y considere el modelo de regresión logística multivariada de STA (variable dependiente) sobre AGE, CAN, CPR, INF y RACE.
La variable RACE está codificada en tres niveles. Prepare una tabla mostrando la codificación de dos variables de diseño necesaria para incluir esta variable en un modelo de regresión logístico. Tome como referencia la categoría White.
Escriba la ecuación general para el modelo de regresión logístico de STA contra AGE, CAN, CPR, INF y RACE y para el logit transformado de este modelo. ¿Cuántos parámetros tiene este modelo?
Escriba una expresión general para la función de verosimilitud y del logaritmo de esta función para el modelo de regresión logístico de (b). ¿Cuántas ecuaciones de verosimilitud hay? Obtenga una expresión para una ecuación de verosimilitud típica´para este problema.
Obtenga las estimaciones de máxima verosimilitud de los parámetros del modelo de regresión logístico de (b). Usando estas estimaciones, escriba las correspondientes ecuaciones para los valores ajustados. Es decir, las probabilidades logísticas estimadas.
Usando los resultados de la salida encontrada en la parte (d), verifique la significancia del coeficiente de AGE. ¿Qué supuestos se necesitan para realizar dicha prueba?
Utilize la estadística de Wald para obtener una aproximación a la significancia individual de los coeficientes de las variables independientes en el modelo. Ajuste un modelo reducido que elimine aquellas variables que no tienen significancia. Ajuste significancia conjunta (condicional) de las variables que quedaron en el modelo.
Usando los resultados de (f), halle un intervalo del \(95\%\) de confianza para todos los coeficientes del modelo. Escriba una interpretación con respecto a los intervalos encontrados para las pendientes.
Obtenga la matriz de covarianzas estimada para el modelo final en (f). Escoja un conjunto de valores para las covariables en ese modelo y calcule el logit y la probabilidad logística estimada para una persona con estas características. Calcule un intervalo del \(95\%\) de confianza para el logit y la probabilidad logística estimada. Interprete la probabilidad estimada y su intervalo de confianza.
9.0.6 Ejercicio 11 a 14
- Los datos PROS corresponden a un estudio realizado pacientes con cáncer de próstata para determinar si las variables medidas en un examen básico pueden ser usadas para predecir si el tumor ha penetrado la cápsula prostática. Los datos fueron recogidos teniendo en cuenta 380 individuos, 153 de los cuales tuvieron un cáncer que penetró la cápsula prostática. En estos datos, entre otras variables, una que se pensó que era particularmente predictiva para la penetración de cápsula es el nivel de antígeno prostático, PSA. Use estos datos y considere el modelo de regresión logística multivariada de CAPSULE (variable dependiente) sobre AGE, RACE, DPROS, DCAPS, PSA, GLEASON y VOL.
La variable DPROS está codificada en cuatro niveles. Prepare una tabla mostrando la codificación de tres variables de diseño necesaria para incluir esta variable en un modelo de regresión logístico. Tome como referencia la categoría No nodule.
La variable DCAPS está codificada en dos niveles (1 y 2). ¿Puede ser usada esta variable en su codificación original o debe ser creada una variable de diseño? Explore esta pregunta comparando las estimaciones de los coeficientes obtenidas en un modelo que contenga DCAPS como se codificó originalmente con aquellas estimaciones encontradas al utilizar una nueva variable de diseño DCAPSnew=DCAPS-1, la cual está codificada con 0 y 1.
Repita las partes (b)-(h) del ejercicio 10.
- Los datos UIS se recogieron con el propósito de comparar dos programas de tratamiento A y B para reducir el abuso de la droga y prevenir sus riesgos. La descipción de los datos se puede ver también aquí. Realize un análisis de regresión logística para encontrar el mejor submodelo logístico, tomando a DFREE como variable dependiente y al resto como variables independientes.
- En los datos CLSLOWBWT, entre otras variables, una que los físicos consideraron importante para el control del peso del bebé (variable dependiente LOW) fue el peso de la madre durante el último periodo menstrual (LWT). Realize un análisis de regresión logística para encontrar el mejor submodelo logístico, tomando a LOW como variable dependiente y al resto como variables independientes.
- Los datos LOWBWT corresponden a un estudio realizado para identificar factores de riesgos asociados a nacimientos de bebés con bajo peso (peso menor que 2.500 gramos). Los datos fueron recogidos teniendo en cuenta 189 mujeres, 59 de las cuales tuvieron bebés con bajo peso y 130 de las cuales tuvieron bebés con peso normal. Realize un análisis de regresión logística para encontrar el mejor submodelo logístico, tomando a LOW como variable dependiente y al resto como variables independientes.
Bibliografía
Consultar el documento RPubs :: Regresión logística (bibliografía).
If you found any ERRORS or have SUGGESTIONS, please report them to my email. Thanks.