Différentes façons d'écrire les termes d'interaction dans lm?

42

J'ai une question sur le meilleur moyen de spécifier une interaction dans un modèle de régression. Considérez les données suivantes:

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")

Deux manières équivalentes de spécifier le modèle avec des interactions sont les suivantes:

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

Ma question est de savoir si je pourrais spécifier l’interaction en considérant une nouvelle variable (rs) avec les mêmes niveaux d’interaction:

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

Quels sont les avantages / inconvénients de cette approche? Et pourquoi les résultats de ces deux approches sont différents?

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
la source
Vous voulez dire que rsc'est défini comme interaction(r, s)?
chl
Peut-être pourriez-vous nous montrer le code qui a créé rsr1s2?
jbowman
Le facteur rs a été défini manuellement (il suffit de coller les facteurs r et s). Voir le jeu de données.
Manuel Ramón
1
Je suppose que c'est lié au chemin dans les variables sont liées voir attr(terms(lm1),"factors")etattr(terms(lm2),"factors")
Galled

Réponses:

8

Les résultats sont différents, car la façon dont lm configure le modèle avec l'interaction diffère de la manière dont il est configuré lorsque vous le configurez vous-même. Si vous regardez le sd résiduel, c'est le même, ce qui indique (de manière non définitive) que les modèles sous-jacents sont les mêmes, ils sont simplement exprimés (aux internes du lm) différemment.

Si vous définissez votre interaction comme étant paste(d$s, d$r)au lieu de paste(d$r, d$s)vos paramètres, les estimations changeront encore, de manière intéressante.

Notez que, dans le résumé de votre modèle pour lm1, l'estimation du coefficient pour ss2 est 4,94 inférieure à celle du résumé pour lm2, avec le coefficient pour rr2: ss2 étant de 4,95 (si vous imprimez avec 3 décimales, la différence disparaîtra). Ceci est une autre indication d'un réarrangement interne des termes.

Je ne vois aucun avantage à le faire vous-même, mais il en existe peut-être un avec des modèles plus complexes dans lesquels vous ne voulez pas un terme d'interaction complet, mais seulement quelques-uns des termes du "croisement" entre deux facteurs ou plus.

Jbowman
la source
Le seul avantage que je vois dans la définition des interactions, comme en lm2, est qu’il est facile d’effectuer des comparaisons multiples pour le terme d’interaction. Ce que je ne comprends pas très bien, c’est pourquoi différents résultats sont obtenus si, en principe, il semble que les deux approches soient les mêmes.
Manuel Ramón
5
X1,X2(1,X1,X2,X1*X2)(X1,X2,X1*X2,(1-X1)*(1-X2)
Par conséquent, bien que différentes, les deux approches sont correctes, n'est-ce pas?
Manuel Ramón
Droite. Mathématiquement, les matrices de variables indépendantes dans les différentes formulations ne sont que des transformations linéaires les unes des autres, de sorte que les estimations de paramètres d'un modèle peuvent être calculées à partir des estimations de paramètres d'un autre si on sait comment les deux modèles ont été réellement configurés.
jbowman
9

Vous pourriez mieux comprendre ce comportement si vous examinez les matrices du modèle.

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

Lorsque vous regardez ces matrices, vous pouvez comparer les constellations de s2=1avec les autres variables (c'est-à-dire quand s2=1, quelles valeurs prennent les autres variables?). Vous verrez que ces constellations diffèrent légèrement, ce qui signifie simplement que la catégorie de base est différente. Tout le reste est essentiellement identique. En particulier, notez que dans votre lm1, le coefficient sur ss2est égal aux coefficients ss2+rsr1s2de lm2, soit 3,82 = 8,76-4,95, sans erreurs d’arrondi.

Par exemple, l'exécution du code suivant vous donne exactement le même résultat que l'utilisation du réglage automatique de R:

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

Cela fournit également une réponse rapide à votre question: la seule raison de changer la configuration des facteurs est de fournir une clarté d'exposition. Prenons l'exemple suivant: supposons que vous régressiez le salaire d'un mannequin pour terminer ses études secondaires en interaction avec un facteur indiquant si vous appartenez à une minorité.

wunege=α+β evous+γ evous*mjenorjety+ε

ββ+γ

coffeinjunky
la source