¿Cómo comparar las curvas de regresión logística?

Question

Tengo una variable dependiente binaria cuya probabilidad depende de una variable independiente continua, es decir, la edad. Las he ajustado en un modelo de regresión logística y tengo los coeficientes de la misma.

El problema es que tengo 3 curvas de este tipo (para 3 grupos diferentes - categóricos) y necesito demostrar que estas tres son estadísticamente diferentes entre sí. ¿Cómo puedo demostrar que una curva de regresión logística no es igual a la otra?

Nota: Los tamaños de las muestras de los tres grupos son diferentes

Answer 1

Yo sugeriría poner los tres conjuntos de datos en uno solo. Incluya en este conjunto de datos una variable indicadora para cada uno de los tres conjuntos de datos. A continuación, ajuste una regresión logística utilizando este conjunto de datos completo, incluyendo un término de interacción de la edad y la variable indicadora. Además, ajuste un modelo de regresión logística que no incluya el término de interacción, sino términos aditivos para la edad y el indicador. Realice una prueba de razón de verosimilitud para evaluar la evidencia de que el término de interacción es significativo. Es decir, compare los modelos

$$ \mathrm{logit}(y)=\beta_0 + \beta_1\cdot \mathrm{age} + \beta_2\cdot \mathrm{g_2} + \beta_3\cdot \mathrm{g_3} + \beta_4\cdot \mathrm{age}\cdot \mathrm{g_2} + \beta_5\cdot \mathrm{age}\cdot \mathrm{g_3} $$

y

$$ \mathrm{logit}(y)=\beta_0 + \beta_1\cdot \mathrm{age} + \beta_2\cdot g_2 + \beta_3\cdot g_3 $$

donde $\mathrm{g_2}$ y $\mathrm{g_3}$ son variables indicadoras (variables ficticias) para los grupos 2 y 3. $g_2$ es $1$ siempre que el grupo sea $2$ y $0$ de lo contrario. $g_3$ es $1$ siempre que el grupo sea $3$ y $0$ otherweise. El grupo 1 sirve como grupo de referencia. Se puede cambiar el grupo de referencia omitiendo la variable ficticia correspondiente en el modelo. En este caso, hemos omitido la variable indicadora del grupo 1 que, por tanto, sirve de grupo de referencia.

La prueba de la razón de verosimilitud pone a prueba las siguientes hipótesis

\begin{align} \mathrm{H}_{0}&: \beta_4 = \beta_5 = 0 \\ \mathrm{H}_{1}&: \text{At least one}\> \beta_j \neq 0, j = 4, 5 \end{align}

Por tanto, comprueba si las diferencias entre la pendiente de la edad del grupo 1 (o el grupo de referencia) y los grupos 2 y 3 difieren significativamente. Sin embargo, no indica si las pendientes de los grupos 2 y 3 difieren (véase la siguiente sección). Si la prueba de razón de verosimilitudes es significativa, tiene pruebas de que al menos dos grupos tienen pendientes diferentes.

Este es un ejemplo en el que se utiliza R . Estoy utilizando un conjunto de datos disponible que modela la probabilidad de admisión en una escuela de posgrado. En lugar de la edad, utiliza el GPA y en lugar del indicador de grupo utiliza los rangos 1-4 que indican el prestigio de la institución. Para tu problema, sólo tienes que cambiar el GPA por la edad y el rango por el indicador de grupo.

dat <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

dat$rank <- factor(dat$rank)

glm_mod <- glm(admit~gpa*rank, family = binomial, data = dat)
glm_mod_noint <- glm(admit~gpa+rank, family = binomial, data = dat)

summary(glm_mod)

# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)  
# (Intercept)  -4.6298     2.4843  -1.864   0.0624 .
# gpa           1.3906     0.7171   1.939   0.0525 .
# rank2         0.9191     2.9618   0.310   0.7563  
# rank3         0.1892     3.2431   0.058   0.9535  
# rank4        -1.1121     4.1319  -0.269   0.7878  
# gpa:rank2    -0.4661     0.8584  -0.543   0.5871  
# gpa:rank3    -0.4567     0.9295  -0.491   0.6232  
# gpa:rank4    -0.1378     1.1999  -0.115   0.9085

summary(glm_mod_noint) # output omitted here

Ahora se ve que la pendiente para el GPA es $1.39$ para las instituciones con rango 1. Para las instituciones con rango 2, la pendiente es $1.391 - 0.466 = 0.925$ la diferencia de pendientes es $-0.466$ . El correspondiente $p$ -valor es $0.587$ , proporcionando poca evidencia de una diferencia en las pendientes. En consecuencia, la diferencia de las pendientes en comparación con las instituciones del rango 1 para el rango 3 es $-0.457$ y $-0.138$ para las instituciones de rango 4. Todos los correspondientes $p$ -son bastante elevados. Hagamos ahora la prueba de la razón de verosimilitud:

# Likelihood ratio test

anova(glm_mod_noint, glm_mod, test = "Chisq")

# Analysis of Deviance Table
# 
# Model 1: admit ~ gpa + rank
# Model 2: admit ~ gpa * rank
#   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
# 1       395     462.88                     
# 2       392     462.49  3  0.38446   0.9434

La alta $p$ -El valor de la prueba de la razón de verosimilitudes indica que hay pocas pruebas de que la influencia del GPA en la probabilidad de admisión difiera entre instituciones de diferentes rangos. Por lo tanto, suponemos que $\beta_5 = \beta_6 = \beta_7 = 0$ en este caso (es decir, los coeficientes de los términos de interacción son todos $0$ ).

Visualicemos nuestros resultados en la escala de probabilidad:

Y ahora en la escala de logaritmos:

Como puede ver, las líneas difieren ligeramente en su pendiente, pero estas diferencias no son estadísticamente significativas. Informalmente, el término de interacción evalúa la evidencia de líneas paralelas. En este caso, podemos suponer que las líneas son paralelas en la escala de logaritmos.

Comparaciones por pares

La prueba de razón de verosimilitud comprueba si todas las diferencias en las pendientes entre el grupo de referencia y los otros grupos son $0$ (véase la hipótesis anterior). Si queremos comprobar todas las diferencias entre pares de pendientes, tenemos que recurrir a una prueba post-hoc. Si hay $n$ grupos, hay $\frac{1}{2}(n - 1)n$ comparaciones por pares. En su caso, hay $n = 3$ grupos, lo que lleva a un total de $3$ comparaciones. En este ejemplo con cuatro grupos (rangos 1-4), hay $6$ comparaciones. Utilizaremos el multcomp para realizar las pruebas:

library(multcomp)

glht_mod <- glht(glm_mod, linfct = c(  "gpa:rank2 = 0"                 # rank 1 vs. rank 2
                                       , "gpa:rank3 = 0"               # rank 1 vs. rank 3
                                       , "gpa:rank4 = 0"               # rank 1 vs. rank 4
                                       , "gpa:rank2 - gpa:rank3 = 0"   # rank 2 vs. rank 3
                                       , "gpa:rank3 - gpa:rank4 = 0"   # rank 3 vs. rank 4
                                       , "gpa:rank2 - gpa:rank4 = 0")) # rank 2 vs. rank 4
summary(dd) # all pairwise tests

# Linear Hypotheses:
#                             Estimate Std. Error z value Pr(>|z|)
# gpa:rank2 == 0             -0.466105   0.858379  -0.543    0.947
# gpa:rank3 == 0             -0.456729   0.929507  -0.491    0.960
# gpa:rank4 == 0             -0.137839   1.199946  -0.115    0.999
# gpa:rank2 - gpa:rank3 == 0 -0.009376   0.756521  -0.012    1.000
# gpa:rank3 - gpa:rank4 == 0 -0.318890   1.129327  -0.282    0.992
# gpa:rank2 - gpa:rank4 == 0 -0.328266   1.071546  -0.306    0.990

Estos $p$ -Los valores se ajustan para las comparaciones múltiples. Hay muy pocas pruebas de cualquier diferencia entre las pendientes. Por último, probemos la hipótesis global de que todo las comparaciones de pendientes por pares son $0$ :

summary(glht_mod, Chisqtest()) # global test

# Global Test:
#    Chisq DF Pr(>Chisq)
# 1 0.3793  3     0.9445

De nuevo, hay muy pocas pruebas de que alguna de las vertientes se diferencie de la otra.