Al usar la función `cbind()` en R para una regresión logística en una tabla de $2 \times 2$, ¿cuál es la forma funcional explícita de la ecuación de regresión?

Question

Supongamos que tengo una tabla de $2 \times 2$ que se ve así:

            Enfermedad       Sin Enfermedad
Tratamiento         55                67
Control           42                34

Me gustaría realizar una regresión logística en R con esta tabla. Entiendo que la forma estándar es usar la función glm con la función cbind en la respuesta. En otras palabras, el código se ve así:

glm(formula = cbind(c(55,67),c(42,34)) ~ as.factor(c(1, 0)), family = binomial())

Me pregunto por qué R nos exige usar la función cbind y por qué simplemente usar proporciones no es suficiente. ¿Hay alguna manera de escribir esto explícitamente como una fórmula? ¿Cómo se vería en la forma de:

$$ log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1X $$

donde $X = 1$ si tenemos tratamiento y $X=0$ para control?

En este momento parece que estoy haciendo una regresión en una matriz para el valor dependiente.

Answer 1

Primero muestro cómo se puede especificar una fórmula utilizando datos agregados con proporciones y pesos. Luego muestro cómo podrías especificar una fórmula después de desagregar tus datos en observaciones individuales.

La documentación en glm indica que:

"Para un GLM binomial se utilizan pesos previos para dar el número de intentos cuando la respuesta es la proporción de éxitos"

Creo nuevas columnas total y proportion_disease en df para 'número de intentos' y 'proporción de éxitos' respectivamente.

library(dplyr)
df <- tibble(treatment_status = c("tratamiento", "sin_tratamiento"),
       enfermedad = c(55, 42),
       no_enfermedad = c(67,34)) %>% 
  mutate(total = no_enfermedad + enfermedad,
         proportion_disease = enfermedad / total) 

modelo_pesado <- glm(proportion_disease ~ treatment_status, 
        data = df, family = binomial("logit"), weights = total)

El enfoque ponderado anterior toma datos agregados y dará la misma solución que el método cbind pero te permite especificar una fórmula. (A continuación es equivalente al método del Póster Original pero cbind(c(55,42), c(67,34)) en lugar de cbind(c(55,67),c(42,34)) para que 'Enfermedad' en lugar de 'Tratamiento' sea la variable de respuesta.)

modelo_cbinded <- glm(cbind(enfermedad, no_enfermedad) ~ 
      treatment_status, data = df, family = binomial("logit"))  

También podrías simplemente desagregar tus datos en observaciones individuales y pasar estos a glm (permitiéndote especificar una fórmula también).

df_expanded <- tibble(estado_enfermedad = c(1, 1, 0, 0), 
           treatment_status = rep(c("tratamiento", "control"), 2)) 
%>%
                        .[c(rep(1, 55), rep(2, 42), rep(3, 67), 
                          rep(4, 34)), ]

modelo_expandido <- glm(estado_enfermedad ~ treatment_status, 
        data = df_expanded, family = binomial("logit"))

Comparemos ahora estos pasando cada modelo a summary. modelo_pesado y modelo_cbinded ambos producen los mismos resultados exactos. modelo_expandido produce los mismos coeficientes y errores estándar, aunque muestra diferentes grados de libertad, deviance, AIC, etc. (correspondiendo con el número de filas/observaciones).

    > lapply(list(modelo_pesado, modelo_cbinded, modelo_expandido), 
             summary)
[[1]]

Llamada:
glm(formula = proportion_disease ~ treatment_status, 
    family = binomial("logit"), 
    data = df, weights = total)

Residuos de deviance: 
[1]  0  0

Coeficientes:
                          Estimado Error Estándar valor z de Pr(>|z|)
(Intercept)                 0.2113     0.2307   0.916    0.360
treatment_statustratamiento  -0.4087     0.2938  -1.391    0.164

(Parámetro de dispersión para la familia binomial tomado como 1)

    Deviance nula: 1.9451e+00  en 1  grados de libertad
Deviance residual: 1.0658e-14  en 0  grados de libertad
AIC: 14.028

Número de iteraciones de Fisher Scoring: 2

[[2]]

Llamada:
glm(formula = cbind(enfermedad, no_enfermedad) ~ treatment_status, 
    family = binomial("logit"), data = df)

Residuos de deviance: 
[1]  0  0

Coeficientes:
                          Estimado Error Estándar valor z de Pr(>|z|)
(Intercept)                 0.2113     0.2307   0.916    0.360
treatment_statustratamiento  -0.4087     0.2938  -1.391    0.164

(Parámetro de dispersión para la familia binomial tomado como 1)

    Deviance nula: 1.9451e+00  en 1  grados de libertad
Deviance residual: 1.0658e-14  en 0  grados de libertad
AIC: 14.028

Número de iteraciones de Fisher Scoring: 2

[[3]]

Llamada:
glm(formula = estado_enfermedad ~ treatment_status, 
     family = binomial("logit"), 
    data = df_expanded)

Residuos de deviance: 
   Min      1Q  Mediana      3Q     Máx  
-1.268  -1.095  -1.095   1.262   1.262  

Coeficientes:
                          Estimado Error Estándar valor z de Pr(>|z|)
(Intercept)                 0.2113     0.2307   0.916    0.360
treatment_statustratamiento  -0.4087     0.2938  -1.391    0.164

(Parámetro de dispersión para la familia binomial tomado como 1)

    Deviance nula: 274.41  en 197  grados de libertad
Deviance residual: 272.46  en 196  grados de libertad
AIC: 276.46

Número de iteraciones de Fisher Scoring: 3

(Ver R bloggers para conversación sobre el parámetro weights en glm en contexto de regresión.)

Answer 2

Distribución binomial versus distribución de Bernoulli

La distribución binomial se puede ver como múltiples distribuciones de Bernoulli con el mismo parámetro de probabilidad $p$. Sus distribuciones de probabilidad difieren hasta una constante.

Para observaciones $x_i \in \lbrace 0,1 \rbrace$, donde $k$ es el número de observaciones que son iguales a $1$

• $1$ Función de masa de probabilidad para el experimento de Bernoulli $$f(x_1) = p^{x_1} (1-p)^{1-x_1}$$
• $n$ Función de masa de probabilidad para $$f(x_1,x_2,\dots,x_n) = \prod p^{x_i} (1-p)^{1-x_i} = p^{k} (1-p)^{n-k}$$ * Distribución binomial de la función de masa de probabilidad $$f(x_1,x_2,\dots,x_n) = f(k,n) = {n \choose k} p^{k} (1-p)^{n-k}$$

**La diferencia es que la distribución binomial kind of 'ignora' el orden de las observaciones $x_i$ y solo mira el número total de $x_i=1$ y $x_i=0$.

Influencia del tamaño de la muestra en la distribución de probabilidad y el error

Por qué simplemente usar proporciones no es suficiente

La estimación de máxima verosimilitud maximiza $p^k(1-p)^{n-k}$ y esto sucede en $p_{max} = k/n$ la proporción de observaciones con $x_i=1$. Entonces, de hecho, para encontrar la máxima verosimilitud la proporción es suficiente.

Sin embargo, cuando deseamos estimar alguna expresión para el error de la estimación (por ejemplo, el error estándar o los intervalos de confianza) entonces el número de observaciones se vuelve importante.

Una vista intuitiva para ver cómo importa el tamaño de la muestra puede ser la siguiente: Veamos la distribución de la muestra del resultado del experimento en función del tamaño de la muestra. En la imagen a continuación trazamos la distribución binomial para $p=0.8$ con diferentes parámetros de tamaño $n$. A la izquierda se ve que las curvas se desplazan hacia la derecha para $n$ más grandes, lo cual es obvio ya que se esperaría más recuentos cuando el tamaño de la muestra $n$ aumenta. Pero la distancia de la tasa esperada $p$ veces $n$ se vuelve relativamente más corta. Se ve esto a la derecha donde dividimos el resultado (los recuentos) por el tamaño de la muestra total $n$, lo que da una tasa relativa. Se ve que la tasa relativa observada es más probable que esté cerca de la tasa real (0.8 en este gráfico) cuando el tamaño $n$ es mayor. (Una pregunta relacionada es: ¿Cómo estimar la probabilidad de que ocurra un evento basado en su recuento?).

Entrada para la función GLM

Me pregunto por qué R nos obliga a usar la función cbind

No es necesario usar la función cbind. Hay 2/3 métodos (dos o tres dependiendo de cómo lo veas: dos formas de presentar los datos y la última se divide en dos formas de alimentarla a la función glm).

• También podemos usar GLM como si estuviera haciendo una variable de Bernoulli. En ese caso, los datos consistirían en $55+67+42+34 = 198$ filas. Una fila para cada recuento individual (cada recuento se puede considerar una variable distribuida de manera Bernoulli)

  Por lo tanto, los datos se verían como una tabla larga

  $$\begin{array}{rrr} \text{ID} & \text{Observación de la enfermedad} & \text{Tratamiento}\\ 1 & 1 & 1\\ 2 & 1 & 1\\ 3 & 1 & 1\\ \vdots & \vdots & \vdots \\ 53 & 1 & 1\\ 54 & 1 & 1\\ 55 & 1 & 1\\ 56 & 0 & 1\\ 57 & 0 & 1\\ 58 & 0 & 1\\ \vdots & \vdots & \vdots \\ 120 & 0 & 1\\ 121 & 0 & 1\\ 122 & 0 & 1\\ 123 & 1 & 0\\ 124 & 1 & 0\\ 125 & 1 & 0\\ \vdots & \vdots & \vdots \\ 162 & 1 & 0\\ 163 & 1 & 0\\ 164 & 1 & 0\\ 165 & 0 & 0\\ 166 & 0 & 0\\ 167 & 0 & 0\\ \vdots & \vdots & \vdots \\ 196 & 0 & 0\\ 197 & 0 & 0\\ 198 & 0 & 0\\ \end{array}$$

   data <- data.frame(ID = 1:198, 
                      observación = 
    c(rep(1,55),rep(0,67),rep(1,42),rep(0,34)),
    tratamiento = c(rep(1,55+67), rep(0,42+34)))
   glm1 <- glm(observation ~ treatment, family = binomial, 
               data = data) 

• El método anterior se usa frecuentemente cuando la variable independiente es continua. En el caso de una variable categórica, podemos agrupar todos los 1s y 0s para una categoría. En este caso, los datos no son 0s y 1s, sino los recuentos de 0s y 1s.

   data <- data.frame(enfermedad = c(55,42),
                      saludable = c(67,34),
                      tratamiento = c(1,0))
  
   ### el método 'cbind' le da a glm los números en cada categoría
   glm2 <- glm(cbind(enfermedad,saludable) ~ tratamiento, 
               family = binomial, data = data)
  
   ### este método le da a glm la fracción y utiliza 'weight' para describir el número de datos
   glm3 <- glm(enfermedad/(enfermedad+saludable) ~ tratamiento, 
     family = binomial, data = data, weight = enfermedad + saludable)

Los resúmenes de las tres llamadas diferentes a glm anteriores proporcionan los mismos resultados. Hay una diferencia en la estadística como la desviación y el AIC, eso se debe a que la función de verosimilitud es diferente por un factor de $n \choose k$ pero eso no es importante porque las medidas relevantes en la comparación de modelos son las tasas de verosimilitud.

> summary(glm1)

Llamada:
glm(fórmula = observation ~ treatment, family = binomial, 
      datos = data)

Residuos de la desviación: 
   Min      1Q  Mediana      3Q     Máx  
-1.268  -1.095  -1.095   1.262   1.262  

Coeficientes:
            Estimado Error estándar Valor z Pr(>|z|)
(Intercepción)   0.2113     0.2307   0.916    0.360
tratamiento    -0.4087     0.2938  -1.391    0.164

(Parámetro de dispersión para familia binomial tomado como 1)

    Desviación nula: 274.41  en 197  grados de libertad
Desviación residual: 272.46  en 196  grados de libertad
AIC: 276.46

Número de iteraciones de 'Fisher Scoring': 3

> summary(glm2)

Llamada:
glm(fórmula = cbind(disease, healthy) ~ treatment, 
    familia = binomial, 
    datos = data)

Residuos de la desviación: 
[1]  0  0

Coeficientes:
            Estimado Error estándar Valor z Pr(>|z|)
(Intercepción)   0.2113     0.2307   0.916    0.360
tratamiento    -0.4087     0.2938  -1.391    0.164

(Parámetro de dispersión para familia binomial tomado como 1)

    Desviación nula: 1.9451e+00  en 1  grados de libertad
Desviación residual: 1.0658e-14  en 0  grados de libertad
AIC: 14.028

Número de iteraciones de 'Fisher Scoring': 2

> summary(glm3)

Llamada:
glm(fórmula = disease/(disease + healthy) ~ treatment, 
    familia = binomial, 
    datos = data, weights = disease + healthy)

Residuos de la desviación: 
[1]  0  0

Coeficientes:
            Estimado Error estándar Valor z Pr(>|z|)
(Intercepción)   0.2113     0.2307   0.916    0.360
tratamiento    -0.4087     0.2938  -1.391    0.164

(Parámetro de dispersión para familia binomial tomado como 1)

    Desviación nula: 1.9451e+00  en 1  grados de libertad
Desviación residual: 1.0658e-14  en 0  grados de libertad
AIC: 14.028

Número de iteraciones de 'Fisher Scoring': 2

>**

Answer 3

Encontré esta discusión muy útil para el análisis que necesito realizar para mi tesis. No estoy segura si lo entendí bien, ¿que los múltiples ensayos de los que hablas son con la misma muestra? En mi caso tengo un grupo de control y un grupo de tratamiento, y cada encuestado responde 4 preguntas donde tiene que elegir entre tren y avión. Ahora quiero analizar la diferencia entre el grupo de control y el grupo de tratamiento en las cuatro opciones combinadas. Así que creo que podría usar el enfoque con un modelo ponderado y la variable dependiente "proporción_tren" (que representa cuántas veces se eligió el tren de las cuatro opciones).

No estoy segura sobre la interpretación de los coeficientes de este modelo. Sé que probablemente son los logaritmos de las razones de probabilidad (o equivalentemente diferencias en los logaritmos de las razones de probabilidad). Pero ¿estas logaritmos de las razones de probabilidad muestran la probabilidad combinada en todos los ensayos (que quiero averiguar), o la probabilidad por ensayo?

Ahora leí en otro lugar que necesito tener en cuenta la correlación dentro de los individuos incluyendo intercepciones aleatorias para las IDs individuales - ya que las cuatro preguntas fueron respondidas por los mismos individuos (pero diferentes individuos en el grupo de control y tratamiento, por supuesto). Ninguno de ustedes aquí incluyó tal intercepción aleatoria para los individuos, así que me pregunto si es necesario o no?