Разные способы написания терминов взаимодействия в лм?

У меня есть вопрос о том, какой способ лучше определить взаимодействие в регрессионной модели. Рассмотрим следующие данные:

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

Два эквивалентных способа указать модель с взаимодействиями:

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

Мой вопрос заключается в том, могу ли я указать взаимодействие с учетом новой переменной (rs) с теми же уровнями взаимодействия:

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

Какие преимущества / недостатки у этого подхода? И почему результаты этих двух подходов отличаются?

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 Мануэль Рамон
источник

Вы имеете в виду, что rsопределяется как interaction(r, s)?

хл

Возможно, вы могли бы показать нам код, который создал rsr1s2?

jbowman

Коэффициент rs был определен вручную (просто вставьте коэффициенты r и s). Смотрите набор данных.

Мануэль Рамон

Я думаю связан путь в переменных relationated Смотри attr(terms(lm1),"factors")иattr(terms(lm2),"factors")

раздражало

Ответы:

Результаты отличаются, потому что способ, которым lm настраивает модель с помощью взаимодействия, отличается от того, как она настраивается, когда вы настраиваете ее самостоятельно. Если вы посмотрите на остаточный sd, это то же самое, что указывает (не окончательно), что лежащие в основе модели одинаковы, просто выражены (внутренним элементам lm) по-разному.

Если вы определите свое взаимодействие, поскольку paste(d$s, d$r)вместо paste(d$r, d$s)ваших параметров оценки снова изменятся интересными способами.

Обратите внимание, что в сводке вашей модели для lm1 оценка коэффициента для ss2 на 4,94 ниже, чем в сводке для lm2, с коэффициентом для rr2: ss2, равным 4,95 (если вы печатаете до 3 десятичных знаков, разница исчезает). Это еще один признак того, что произошла внутренняя перестановка терминов.

Я не могу думать о каком-либо преимуществе делать это самостоятельно, но может быть одна с более сложными моделями, где вам не нужен термин полного взаимодействия, а вместо этого - только некоторые из терминов в «кресте» между двумя или более факторами.

jbowman
источник

Единственное преимущество, которое я вижу для определения взаимодействий, как в lm2, состоит в том, что легко выполнить несколько сравнений для термина взаимодействия. Что я не совсем понимаю, так это почему разные результаты получаются, если в принципе кажется, что два подхода одинаковы.

Мануэль Рамон

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

Поэтому, хотя они разные, оба подхода верны, не так ли?

Мануэль Рамон

Правильно. Математически матрицы независимых переменных в различных формулировках являются просто линейными преобразованиями друг друга, поэтому оценки параметров одной модели могут быть вычислены из оценок параметров другой, если известно, как на самом деле были установлены две модели.

jbowman

Вы можете лучше понять это поведение, если посмотрите на матрицы моделей.

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

Когда вы смотрите на эти матрицы, вы можете сравнить созвездия s2=1с другими переменными (то есть когда s2=1, какие значения принимают другие переменные?). Вы увидите, что эти созвездия немного отличаются, что означает, что базовая категория отличается. Все остальное по сути то же самое. В частности, обратите внимание, что по вашему lm1, коэффициент наss2 равен коэффициенты ss2+rsr1s2из lm2, т.е. 3,82 = 8.76-4.95, короткие ошибок округления.

Например, выполнение следующего кода даст вам точно такой же вывод, как при использовании автоматической установки R:

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

Это также дает быстрый ответ на ваш вопрос: единственная причина, по которой нужно изменить способ установки факторов, заключается в обеспечении ясности изложения. Рассмотрим следующий пример. Предположим, что вы регрессируете заработную плату на манекене для окончания средней школы, взаимодействуя с фактором, указывающим, принадлежите ли вы к меньшинству.

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

$\beta$ $\beta+\gamma$

coffeinjunky
источник