Różne sposoby pisania warunków interakcji w lm?

Mam pytanie, w jaki sposób najlepiej określić interakcję w modelu regresji. Rozważ następujące dane:

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

Dwa równoważne sposoby określania modelu z interakcjami to:

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

Moje pytanie brzmi: czy mógłbym określić interakcję, biorąc pod uwagę nową zmienną (rs) o tych samych poziomach interakcji:

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

Jakie zalety / wady mają to podejście? I dlaczego wyniki tych dwóch podejść są różne?

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

r regression interaction

— Manuel Ramón
źródło

Masz na myśli rsto, że jest zdefiniowane jako interaction(r, s)?

— chl.

Być może mógłbyś nam pokazać kod, który stworzył rsr1s2?

— jbowman,

Współczynnik rs został zdefiniowany ręcznie (wystarczy wkleić współczynniki ris). Zobacz zestaw danych.

— Manuel Ramón

attr(terms(lm1),"factors")attr(terms(lm2),"factors")

— Wydaje

Odpowiedzi:

Wyniki są różne, ponieważ sposób, w jaki lm konfiguruje model z interakcją, różni się od tego, jak jest konfigurowany podczas samodzielnej konfiguracji. Jeśli spojrzysz na resztkowe sd, jest to to samo, co wskazuje (nie definitywnie), że podstawowe modele są takie same, po prostu inaczej wyrażone (do wewnętrznych elementów lm).

Jeśli zdefiniujesz interakcję jako paste(d$s, d$r)zamiast paste(d$r, d$s)parametrów, szacunki zmienią się ponownie, w ciekawy sposób.

Zwróć uwagę, że w podsumowaniu modelu dla lm1 oszacowanie współczynnika dla ss2 jest o 4,94 niższe niż w podsumowaniu dla lm2, przy współczynniku dla rr2: ss2 równym 4,95 (jeśli drukujesz z dokładnością do 3 miejsc po przecinku, różnica zniknie). To kolejna wskazówka, że nastąpiła wewnętrzna zmiana warunków.

Nie mogę wymyślić żadnej korzyści z robienia tego samemu, ale może istnieć jeden z bardziej złożonymi modelami, w którym nie chcesz pełnego terminu interakcji, ale zamiast tego tylko niektóre z terminów „krzyżujących” dwa lub więcej czynników.

— łucznik
źródło

Jedyną zaletą, jaką widzę, aby zdefiniować interakcje, tak jak w lm2, jest to, że łatwo jest wykonać wiele porównań dla terminu interakcji. Nie do końca rozumiem, dlaczego uzyskuje się różne wyniki, jeśli w zasadzie wydaje się, że oba podejścia są takie same.

— Manuel Ramón

x_{1}, x_{2}

$x_1, x_2$

(1, x_{1}, x_{2}, x_{1} * x_{2})

$(1, x_1, x_2, x_1*x_2)$

(x_{1}, x_{2}, x_{1} * x_{2}, (1 - x_{1}) * (1 - x_{2})

$(x_1, x_2, x_1*x_2, (1-x_1)*(1-x_2)$

Dlatego, choć różne, oba podejścia są poprawne, prawda?

— Manuel Ramón

Dobrze. Matematycznie macierze zmiennych niezależnych w różnych formułach są po prostu liniowymi transformacjami względem siebie, więc oszacowania parametrów jednego modelu można obliczyć na podstawie oszacowań parametrów innego modelu, jeśli wiadomo, w jaki sposób dwa modele zostały faktycznie ustawione.

— jbowman,

Możesz lepiej zrozumieć to zachowanie, jeśli spojrzysz na matryce modelu.

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

Kiedy patrzysz na te macierze, możesz porównać konstelacje s2=1z innymi zmiennymi (tj. Kiedy s2=1, jakie wartości przyjmują inne zmienne?). Zobaczysz, że te konstelacje różnią się nieznacznie, co oznacza po prostu, że kategoria podstawowa jest inna. Cała reszta jest zasadniczo taka sama. W szczególności należy zwrócić uwagę, że w lm1współczynnik na ss2równe współczynniki ss2+rsr1s2o lm2, tj 3,82 = 8.76-4.95, krótkie zaokrąglania błędów.

Na przykład wykonanie następującego kodu daje dokładnie takie same wyniki, jak przy użyciu automatycznego ustawienia R:

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

Zapewnia to również szybką odpowiedź na twoje pytanie: naprawdę jedynym powodem zmiany sposobu konfiguracji czynników jest zapewnienie jasności ekspozycyjnej. Rozważ następujący przykład: Załóżmy, że regresujesz wynagrodzenie dla manekina, aby ukończyć szkołę średnią, mając do czynienia z czynnikiem wskazującym, czy należysz do mniejszości.

$wage = \alpha + \beta \ edu + \gamma \ edu*minority + \epsilon$

$\beta$ $\beta+\gamma$

— coffeinjunky
źródło