¿Por qué algunas estimaciones de regresión difieren por un cambio de signo, pero otras no, cuando cambio el nivel de referencia?

Question

Supongamos que tengo un resultado continuo y y dos predictores factoriales, cada uno con dos niveles. Uno de mis predictores categóricos, drug puede tener dos niveles ("A" o "B"), el otro es smokeYes . Cuando ejecuto un modelo de regresión puedo elegir el nivel de base o de referencia de drug para ser "A", como hice en model1 :

set.seed(123)
y<-rnorm(100, 100, 10)
drug.ab<-factor(sample(c("A", "B"), 100, T), levels=c("A", "B"))
drug.ba<-factor(drug.ab, levels=c("B", "A"))
smoke<-factor(sample(c("Yes", "No"), 100, T), levels=c("No", "Yes"))

#model1:
coef(summary(lm(y~drug.ab*smoke)))
                     Estimate Std. Error    t value     Pr(>|t|)
(Intercept)       100.7484158   2.065091 48.7864379 1.465848e-69
drug.abB            0.9030541   2.796146  0.3229639 7.474250e-01
smokeYes           -0.8693598   2.632484 -0.3302431 7.419359e-01
drug.abB:smokeYes   0.8709116   3.746684  0.2324487 8.166844e-01

O puedo establecer la línea de base en "B", como hice en model2 :

#model2:
coef(summary(lm(y~drug.ba*smoke)))
                       Estimate Std. Error       t value     Pr(>|t|)
(Intercept)       101.651469922   1.885161 53.9218978856 1.377147e-73
drug.baA           -0.903054145   2.796146 -0.3229638818 7.474250e-01
smokeYes            0.001551843   2.666021  0.0005820821 9.995368e-01
drug.baA:smokeYes  -0.870911601   3.746684 -0.2324486531 8.166844e-01

Mi pregunta es por qué la estimación de smokeYes difiere entre model1 et model2 ? ¿Por qué no se diferencia por un cambio de signo como drug.baA ¿y el término de interacción?

Answer 1

Permíteme preparar un ejemplo sencillo para explicarte el concepto, y luego podemos cotejarlo con tus coeficientes.

Tenga en cuenta que al incluir tanto la variable ficticia "A/B" como el término de interacción, está dando a su modelo la flexibilidad de ajustar un intercepto (utilizando la variable ficticia) y una pendiente (utilizando la interacción) diferentes en los datos "A" y en los datos "B". En lo que sigue, realmente no importa si el otro predictor $x$ es una variable continua o, como en su caso, otra variable ficticia. Si hablo en términos de su "intercepción" y "pendiente", esto puede interpretarse como "nivel cuando la variable ficticia es cero" y "cambio de nivel cuando la variable ficticia pasa de $0$ a $1$ " si lo prefiere.

Supongamos que el modelo ajustado por OLS en los datos "A" solamente es $\hat y = 12 + 5x$ y en los datos "B" solo es $\hat y = 11 + 7x$ . Los datos podrían tener este aspecto:

Supongamos ahora que tomamos "A" como nivel de referencia y utilizamos una variable ficticia $b$ para que $b=1$ para las observaciones del Grupo B, pero $b=0$ en el Grupo A. El modelo ajustado en todo el conjunto de datos es

$$\hat y_i = \hat \beta_0 + \hat \beta_1 x_i + \hat \beta_2 b_i + \hat \beta_3 x_ib_i$$

Para las observaciones del grupo A tenemos $\hat y_i = \hat \beta_0 + \hat \beta_1 x_i$ y podemos minimizar su suma de residuos al cuadrado estableciendo $\hat \beta_0 = 12$ et $\hat \beta_1 = 5$ . Para los datos del Grupo B, $\hat y_i = (\hat \beta_0 + \hat \beta_2) + (\hat \beta_1 + \hat \beta_3) x_i$ y podemos minimizar su suma de residuos al cuadrado tomando $\hat \beta_0 + \hat \beta_2 = 11$ et $\hat \beta_1 + \hat \beta_3 = 7$ . Está claro que podemos minimizar la suma de los residuos al cuadrado en la regresión global minimizando las sumas para ambos grupos, y que esto se puede conseguir estableciendo $\hat \beta_0 = 12$ et $\hat \beta_1 = 5$ (del Grupo A) y $\hat \beta_2 = -1$ et $\hat \beta_3 = 2$ (ya que los datos "B" deberían tener un intercepto uno inferior y una pendiente dos superior). Obsérvese cómo la presencia de un término de interacción era necesaria para que tuviéramos la flexibilidad suficiente para minimizar la suma de los residuos al cuadrado para ambos grupos a la vez . Mi modelo ajustado será:

$$\hat y_i = 12 + 5 x_i - 1 b_i +2 x_i b_i$$

Cambia todo esto para que "B" sea el nivel de referencia y $a$ es una variable ficticia que codifica el Grupo A. ¿Puede ver que ahora debo ajustar el modelo

$$\hat y_i = 11 + 7 x_i + 1 a_i -2 x_i a_i$$

Es decir, tomo el intercepto ( $11$ ) y la pendiente ( $7$ ) de mi grupo "B" de referencia, y utilizar el término ficticio y de interacción para ajustarlos para mi grupo "A". Estos ajustes esta vez son en la dirección inversa (necesito un intercepto más alto y una pendiente dos baja ) por lo que los signos se invierten en comparación con cuando tomé "A" como grupo de referencia, pero debería estar claro por qué los otros coeficientes no han cambiado simplemente de signo.

---

Comparemos eso con su resultado. En una notación similar a la anterior, su primer modelo ajustado con la línea de base "A" es:

$$\hat y_i = 100.7484158 + 0.9030541 b_i -0.8693598 x_i + 0.8709116 x_i b_i$$

Su segundo modelo ajustado con la línea de base "B" es:

$$\hat y_i = 101.651469922 -0.903054145 a_i + 0.001551843 x_i -0.870911601 x_i a_i$$

En primer lugar, comprobemos que estos dos modelos van a dar los mismos resultados. Pongamos $b_i = 1 - a_i$ en la primera ecuación, y obtenemos

$$\hat y_i = 100.7484158 + 0.9030541 (1-a_i) -0.8693598 x_i + 0.8709116 x_i (1-a_i)$$

Esto se simplifica a:

$$\hat y_i = (100.7484158 + 0.9030541) - 0.9030541 a_i + (-0.8693598 + 0.8709116) x_i - 0.8709116 x_i a_i$$

Una rápida operación aritmética confirma que este modelo es el mismo que el segundo modelo ajustado; además, ahora debería quedar claro qué coeficientes han cambiado de signo y qué coeficientes se han ajustado simplemente a la otra línea de base.

En segundo lugar, veamos cuáles son los diferentes modelos ajustados en los grupos "A" y "B". Su primer modelo da inmediatamente $\hat y_i = 100.7484158 -0.8693598 x_i$ para el grupo "A", y su segundo modelo da inmediatamente $\hat y_i = 101.651469922 + 0.001551843 x_i$ para el grupo "B". Puede comprobar que el primer modelo da el resultado correcto para el grupo "B" sustituyendo $b_i = 1$ en su ecuación; el álgebra, por supuesto, funciona de la misma manera que el ejemplo más general anterior. Del mismo modo, se puede comprobar que el segundo modelo da el resultado correcto para el grupo "A" estableciendo $a_i = 1$ .

En tercer lugar, como en su caso el otro regresor era también una variable ficticia, le sugiero que calcule las medias condicionales ajustadas para las cuatro categorías ("A" con $x=0$ , "A" con $x=1$ , "B" con $x=0$ , "B" con $x=1$ ) en ambos modelos y compruebe que entiende por qué coinciden. Estrictamente hablando, esto es innecesario, ya que hemos realizado el álgebra más general anterior para mostrar que los resultados serán consistentes incluso si $x$ es continua pero creo que sigue siendo un ejercicio valioso. No voy a entrar en detalles, ya que la aritmética es sencilla y está más en línea con el espíritu de la muy agradable respuesta de JonB. Un punto clave que hay que entender es que, sea cual sea el grupo de referencia que se utilice, el modelo tiene suficiente flexibilidad para ajustar cada media condicional por separado. (Aquí es donde marca la diferencia que su $x$ es una variable ficticia para un factor binario en lugar de una variable continua - con predictores continuos generalmente no esperamos que la media condicional estimada $\hat y$ para que coincida con la media de la muestra para cada combinación observada de predictores). Calcule la media de la muestra para cada una de esas cuatro combinaciones de categorías, y debería encontrar que coinciden con sus medias condicionales ajustadas.

Código R para dibujar y explorar los modelos ajustados, las predicciones $\hat y$ y medios de grupo

#Make data set with desired conditional means
data.df <- data.frame(
  x = c(0,0,0,        1,1,1,        0,0,0,        1,1,1),
  b = c(0,0,0,        0,0,0,        1,1,1,        1,1,1),
  y = c(11.8,12,12.2, 16.8,17,17.2, 10.8,11,11.2, 17.8,18,18.2)
)
data.df$a <- 1 - data.df$b

baselineA.lm <- lm(y ~ x * b, data.df)
summary(baselineA.lm) #check this matches y = 12 + 5x - 1b + 2xb

baselineB.lm <- lm(y ~ x * a, data.df)
summary(baselineB.lm) #check this matches y = 11 + 7x + 1a - 2xa

fitted(baselineA.lm)
fitted(baselineB.lm) #check the two models give the same fitted values for y...
with(data.df, tapply(y, interaction(x, b), mean)) #...which are the group sample means

colorSet <- c("red", "blue")
symbolSet <- c(19,17)
with(data.df, plot(x, y, yaxt="n", col=colorSet[b+1], pch=symbolSet[b+1],
                   main="Response y against other predictor x",
                   panel.first = {
                     axis(2, at=10:20)
                     abline(h = 10:20, col="gray70")
                     abline(v = 0:1,  col="gray70")
                   }))
abline(lm(y ~ x, data.df[data.df$b==0,]), col=colorSet[1])
abline(lm(y ~ x, data.df[data.df$b==1,]), col=colorSet[2])
legend(0.1, 17, c("Group A", "Group B"), col = colorSet,
       pch = symbolSet, bg = "gray95")

Answer 2

Eso tiene que ver con la forma en que se define la intercepción. En el primer ejemplo, el intercepto se define como aquellos que no fuman y que tienen la droga A. Los fumadores, que también tienen la droga A, tendrán un valor de 100,75 - 0,87 = 99,9 mientras que los fumadores que tienen la droga B tendrán un valor de 100,75 + 0,90 - 0,87 + 0,87 = 101,65.

En el segundo ejemplo, el intercepto se define como aquellos que no fuman y tienen la droga B. Los fumadores con la droga B tendrán entonces un valor de 101,65 + 0,001 = 101,65, y los fumadores con la droga A tendrán un valor de 100,65 - 0,90 + 0,001-0,87 = 99,9.

Así que todo se suma, es sólo una cuestión de cómo se define el intercepto, es decir, el nivel cuando todos los factores se establecen en la categoría de referencia.