¿Diferentes formas de escribir términos de interacción en lm?

42

Tengo una pregunta sobre cuál es la mejor manera de especificar una interacción en un modelo de regresión. Considere los siguientes datos:

d <- structure(list(r = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
     1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("r1","r2"),
     class = "factor"), s = structure(c(1L, 1L, 1L, 1L, 1L, 
     2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), 
    .Label = c("s1","s2"), class = "factor"), rs = structure(c(1L, 1L,
     1L,1L, 1L,2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L),
    .Label = c("r1s1","r1s2", "r2s1", "r2s2"), class = "factor"), 
     y = c(19.3788027518437, 23.832287726332, 26.2533235300492,
     15.962906892112, 24.2873740664331, 28.5181676764727, 25.2757801195961,
     25.3601044326474, 25.3066440027202, 24.3298865128677, 32.5684219007394,
     31.0048406654209, 31.671238316086, 34.1933764518288, 36.8784821769123,
     41.6691435168277, 40.4669714825801, 39.2664137501106, 39.4884849591932,
     49.247505535468)), .Names = c("r","s", "rs", "y"), 
     row.names = c(NA, -20L), class = "data.frame")

Dos formas equivalentes de especificar el modelo con interacciones son:

lm0 <- lm(y ~ r*s, data=d)
lm1 <- lm(y ~ r + s + r:s, data=d)

Mi pregunta es si podría especificar la interacción considerando una nueva variable (rs) con los mismos niveles de interacción:

lm2 <- lm(y ~ r + s + rs, data=d)

¿Qué ventajas / desventajas tiene este enfoque? ¿Y por qué los resultados de estos dos enfoques son diferentes?

summary(lm1)

lm(formula = y ~ r + s + r:s, data = d, x = TRUE)
            coef.est coef.se
(Intercept) 21.94     1.46  
rr2         11.32     2.07  
ss2          3.82     2.07  
rr2:ss2      4.95     2.92  
---
n = 20, k = 4
residual sd = 3.27, R-Squared = 0.87


summary(lm2)

lm(formula = y ~ r + s + rs, data = d, x = TRUE)
            coef.est coef.se
(Intercept) 21.94     1.46  
rr2         11.32     2.07  
ss2          8.76     2.07   # ss2 coef is different from lm1
rsr1s2      -4.95     2.92  
---
n = 20, k = 4
residual sd = 3.27, R-Squared = 0.87
Manuel Ramón
fuente
¿Quieres decir que rsse define como interaction(r, s)?
chl
¿Quizás podría mostrarnos el código que creó rsr1s2?
jbowman
El factor rs se definió manualmente (simplemente pegue los factores r y s). Ver el conjunto de datos.
Manuel Ramón
1
Supongo que está relacionado con la forma en que las variables están relacionadas ver attr(terms(lm1),"factors")yattr(terms(lm2),"factors")
Galled 02 de

Respuestas:

8

Los resultados son diferentes porque la forma en que lm configura el modelo con la interacción es diferente de cómo se configura cuando lo configura usted mismo. Si observa el sd residual, es el mismo, lo que indica (no definitivamente) que los modelos subyacentes son los mismos, solo expresados ​​(a las partes internas de lm) de manera diferente.

Si define su interacción como en paste(d$s, d$r)lugar de paste(d$r, d$s)sus parámetros, las estimaciones cambiarán nuevamente, de maneras interesantes.

Observe cómo en el resumen de su modelo para lm1 el coeficiente estimado para ss2 es 4.94 más bajo que en el resumen para lm2, con el coeficiente para rr2: ss2 es 4.95 (si imprime con 3 decimales, la diferencia desaparece). Esta es otra indicación de que se ha producido una reorganización interna de los términos.

No se me ocurre ninguna ventaja para hacerlo usted mismo, pero puede haber uno con modelos más complejos en los que no desea un término de interacción completo, sino solo algunos de los términos en el "cruce" entre dos o más factores.

jbowman
fuente
La única ventaja que veo para definir las interacciones como en lm2 es que es fácil realizar múltiples comparaciones para el término de interacción. Lo que no entiendo es por qué se obtienen resultados diferentes si, en principio, parece que los 2 enfoques son los mismos.
Manuel Ramón
55
Los enfoques son los mismos, pero las parametrizaciones exactas del modelo que se estima son diferentes, por lo que los resultados parecen diferentes. Considere un modelo con dos regresores binarios y una interacción. Tiene cuatro categorías, pero puede escribir el modelo de varias maneras diferentes, por ejemplo, que 1 sea un término constante, con variables o u otros. Las variables son solo combinaciones lineales entre sí. Las estimaciones de los coeficientes serán diferentes, pero el modelo es realmente el mismo.x1,x2(1,x1,x2,x1x2)(x1,x2,x1x2,(1x1)(1x2)
jbowman
Por lo tanto, aunque diferentes, ambos enfoques son correctos, ¿no?
Manuel Ramón
Correcto. Matemáticamente, las matrices de variables independientes en las diversas formulaciones son solo transformaciones lineales entre sí, por lo que las estimaciones de parámetros de un modelo se pueden calcular a partir de las estimaciones de parámetros de otro si se sabe cómo se configuraron realmente los dos modelos.
jbowman
9

Es posible que comprenda mejor este comportamiento si observa las matrices modelo.

 model.matrix(lm1 <- lm(y ~ r*s, data=d))
 model.matrix(lm2 <- lm(y ~ r + s + rs, data=d))

Cuando observa estas matrices, puede comparar las constelaciones de s2=1con las otras variables (es decir s2=1, ¿ cuándo , qué valores toman las otras variables?). Verá que estas constelaciones difieren ligeramente, lo que significa que la categoría base es diferente. Todo lo demás es esencialmente lo mismo. En particular, tenga en cuenta que en su lm1, el coeficiente en ss2es igual a los coeficientes ss2+rsr1s2de lm2, es decir, 3.82 = 8.76-4.95, menos errores de redondeo.

Por ejemplo, ejecutar el siguiente código le proporciona exactamente el mismo resultado que usar la configuración automática de R:

  d$rs <- relevel(d$rs, "r1s1")
  summary(lm1 <- lm(y~ factor(r) + factor(s) + factor(rs), data=d))

Esto también proporciona una respuesta rápida a su pregunta: la única razón para cambiar la forma en que se configuran los factores es proporcionar claridad expositiva. Considere el siguiente ejemplo: suponga que regresa el salario de un muñeco para completar la escuela secundaria y que interactúa con un factor que indica si pertenece a una minoría.

Es decir:wage=α+β edu+γ eduminority+ϵ

Si dicho factor minoritario toma el valor 1 si pertenece a una minoría, el coeficiente puede interpretarse como una diferencia salarial para las personas no minoritarias que han completado la escuela secundaria. Si este es su coeficiente de interés, entonces debe codificarlo como tal. De lo contrario, suponga que el factor minoritario toma el valor de 1 si no pertenece a una minoría. Luego, para ver cuánto más ganan las personas no pertenecientes a minorías cuando terminan la escuela secundaria, tendrías que calcular "manualmente" . Sin embargo, tenga en cuenta que toda la información está contenida en las estimaciones y que los resultados sustanciales no cambian al configurar los factores de manera diferente.β + γββ+γ

coffeinjunky
fuente