cómo interpretar la dependencia de clases que no son la clase de referencia en un modelo lineal

Question

Si ejecutamos los tres códigos siguientes:

n <- 250
df <- rbind(
  data.frame(cat=1, x=rnorm(n, 0), y=rnorm(n, 1)),
  data.frame(cat=2, x=rnorm(n, 0), y=rnorm(n, 0)),
  data.frame(cat=3, x=rnorm(n, 0), y=rnorm(n, 0))
)
df$cat <- as.factor(df$cat)
df$cat <- relevel(df$cat, ref = "1")

summary(lm(y ~ cat, df))

Aquí dirá que intercept, cat2 y cat3 son estadísticamente significativos al intentar predecir y. Pero no indica que cat=2 y cat=3 sean realmente redundantes.

n <- 250
df <- rbind(
  data.frame(cat=1, x=rnorm(n, 0), y=rnorm(n, 1)),
  data.frame(cat=2, x=rnorm(n, 0), y=rnorm(n, 0)),
  data.frame(cat=3, x=rnorm(n, 0), y=rnorm(n, 0))
)
df$cat <- as.factor(df$cat)
df$cat <- relevel(df$cat, ref = "2")

summary(lm(y ~ cat, df))

Aquí sólo cat1 es significativo.

n <- 250
df <- rbind(
  data.frame(cat=1, x=rnorm(n, 0), y=rnorm(n, 1)),
  data.frame(cat=2, x=rnorm(n, 0), y=rnorm(n, 0)),
  data.frame(cat=3, x=rnorm(n, 0), y=rnorm(n, 0))
)
df$cat <- as.factor(df$cat)
df$cat <- relevel(df$cat, ref = "3")

summary(lm(y ~ cat, df))

Lo mismo digo (lo cual tiene sentido).

¿Tenemos que ejecutar varias veces el modelo lineal para ver que cat2 y cat3 no necesitan utilizarse ambas para predecir y? ¿Y si sólo el intercepto es significativo, qué significa eso? No soy capaz de reproducir un caso en el que eso suceda y ocurre en mi conjunto de datos. ¿Significa eso que no necesitamos ninguna de las 3 variables categóricas para predecir y? ¿Por qué no serían todos los valores p no significativos en ese caso?

Answer 1

Se trata de una confusión frecuente derivada de la forma en que suelen mostrarse las estimaciones de coeficientes y los valores p para los predictores categóricos multinivel. Con esta codificación de tratamiento por defecto, el intercepto es el resultado en el nivel de referencia y los coeficientes son el diferencias de ese resultado asociado a los otros niveles.

La elección del nivel de referencia importa en términos de visualización, pero no para el modelo fundamental. Si cambia el nivel de referencia, tanto el intercepto como las estimaciones de los demás coeficientes cambiarán .

En este caso concreto, puede obtener una visualización más clara si omite la intercepción. (Fijé la semilla en 2552 antes de continuar con su código para generar df, y mantuve el nivel de referencia en "1").

lm2 <- lm(y ~ -1 + cat, df)

En lugar de centrarse en los coeficientes individuales y los valores p, examine la importancia global del predictor categórico. Con estos datos, la norma anova() es una buena elección (sobre todo porque se trata de un ANOVA equilibrado en lm forma). El predictor es altamente significativo, como se pretendía.

anova(lm2)
# Analysis of Variance Table
# 
# Response: y
#            Df Sum Sq Mean Sq F value    Pr(>F)    
# cat         3 269.47  89.824  86.183 < 2.2e-16
# Residuals 747 778.56   1.042                      

Una vez omitido el intercepto, cada uno de los niveles se muestra con su propio valor medio y la significación de su diferencia respecto a 0.

summary(lm2)

## some output lines omitted

# Coefficients:
#      Estimate Std. Error t value Pr(>|t|)    
# cat1  1.03630    0.06457  16.050   <2e-16
# cat2 -0.05519    0.06457  -0.855    0.393    
# cat3  0.03040    0.06457   0.471    0.638    

Sin embargo, este truco sólo funciona con un único predictor categórico.

En general, con la codificación del tratamiento, el intercepto representa una situación en la que todos los predictores categóricos están en sus niveles de referencia y los predictores continuos están en 0. Por tanto, la "significación" del intercepto es simplemente si el resultado es significativamente diferente de 0 para esa elección concreta de niveles de referencia y centrado de predictores continuos . En muchas circunstancias, esa "importancia" no es muy útil por sí sola.

Del mismo modo, la "significación" de un coeficiente para un nivel no de referencia es si su diferencia del resultado en el nivel de referencia puede distinguirse de 0. Así que eso también depende de la elección del nivel de referencia.

No obstante, las predicciones de resultados para cualquier nivel de un predictor categórico serán las mismas independientemente del nivel de referencia o de las opciones de centrado del predictor. A menudo es aconsejable ignorar los valores p de los coeficientes en esta situación y fijarse directamente en las predicciones o en las diferencias entre predicciones.

Les emmeans paquete proporciona una forma general de mostrar resúmenes de modelos más útiles de ese tipo. Funciona con el modelo original para mostrar estimaciones particulares o comparaciones de interés.

library(emmeans)
emmeans(lm2,pairwise~cat)
# $emmeans
#  cat  emmean     SE  df lower.CL upper.CL
#  1    1.0363 0.0646 747   0.9095   1.1631
#  2   -0.0552 0.0646 747  -0.1819   0.0716
#  3    0.0304 0.0646 747  -0.0964   0.1572
# 
# Confidence level used: 0.95 
# 
# $contrasts
#  contrast estimate     SE  df t.ratio p.value
#  1 - 2      1.0915 0.0913 747  11.953  <.0001
#  1 - 3      1.0059 0.0913 747  11.016  <.0001
#  2 - 3     -0.0856 0.0913 747  -0.937  0.6168
# 
# P value adjustment: tukey method for comparing a family of 3 estimates 
# 

Les $emmeans son las estimaciones de los resultados para cada nivel, el $contrasts son las estimaciones de las diferencias entre cada par de niveles de predictor, con una corrección del valor p para comparaciones múltiples. Obtendría los mismos resultados si utilizara modelos de los mismos datos incluyendo el intercepto, o con cualquier otra elección de nivel de referencia.