Определите, какой главный эффект относится к какому взаимодействию

Предположим, у нас есть следующая модель:

set.seed(1)
d <- data.frame(a = gl(4, 1, 64), a4 = sample(4, 64, TRUE), 
                x = rnorm(64), y = rnorm(64))

l <- lm(y ~ a4 + a * x, d)

За взаимодействие x:a получу 3 коэффициента x:a2, x:a3, x:a4. Теперь я хочу определить, какие коэффициенты являются соответствующими основными эффектами, связанными с этим взаимодействием, то есть x, a2, a3 и a4.

Моя идея состояла в том, чтобы использовать strsplit для взаимодействий и получить соответствующие основные эффекты:

(atoms <- strsplit(names(coef(l))[7:9], ":"))
# [[1]]
# [1] "a2" "x" 
# [[2]]
# [1] "a3" "x" 
# [[3]]
# [1] "a4" "x"

Все идет нормально. Но теперь я хотел бы получить значение соответствующего основного эффекта. Хотя это прямолинейно для x, a2, a3 (поскольку это уникальные имена), я изо всех сил пытаюсь понять, как я могу сделать это с a4:

lapply(atoms, function(.) coef(l)[.])
# [[1]]
#        a2         x 
# 0.3630732 0.2136751 
# [[2]]
#         a3          x 
# 0.04153299 0.21367510 
# [[3]]
#         a4          x 
# 0.04765737 0.21367510

Результат для a4 неверен, потому что это основной эффект, связанный с переменной a4, а не с фиктивным закодированным фактором 4 фактора a.

Итак, модель, которую я показывал, является допустимой моделью в R, но имена коэффициентов неоднозначны. Итак, есть ли другой способ сделать правильное сопоставление между коэффициентами взаимодействия и соответствующими основными эффектами?


r interaction
person thothal    schedule 23.03.2016    source источник
comment
Почему бы просто не переименовать столбец a4 заранее? names(d)[which(names(d)=="a4")] <- "a4_".   -  person lukeA    schedule 23.03.2016
comment
Или взять последний a4 из всех a4: lapply(atoms, function(x) coef(l)[!duplicated(names(coef(l)), fromLast = T)][x]) .   -  person lukeA    schedule 23.03.2016
comment
Что ж, на самом деле это не решает проблему устойчивым образом, так как зависит от порядка параметров в модели (которые я не могу контролировать). Я мог бы переделать модель, чтобы обеспечить уникальные имена, но мне это кажется излишним. Единственный способ установить связь через названия эффектов?   -  person thothal    schedule 23.03.2016

Ответы (1)


arrow_upward
1
arrow_downward

Вы можете использовать компонент assign объекта lm:

l$assign
[1] 0 1 2 2 2 3 4 4 4

Это сопоставляет коэффициенты с расширенной формулой a4 + a + x + a:x.

См. help("model.matrix") документацию по компоненту assign.

Изменить:

Чтобы расширить мой ответ, вы можете сделать это:

terms <- labels(terms(l))
coef(l)[l$assign %in% which(terms %in% strsplit("a:x", ":", fixed = TRUE)[[1]])]

#        a2         a3         a4          x 
#0.36307321 0.04153299 0.23383125 0.21367510
person Roland    schedule 25.03.2016
comment
Да, я знаю этот слот, но как он помогает определить соответствующие основные эффекты? Я могу использовать assign, чтобы определить, что коэффициенты 3:5 принадлежат a, но как определить соответствующие основные эффекты? - person thothal; 29.03.2016
comment
УХ ТЫ! Вы сделали это - большое спасибо, именно то, что я искал! - person thothal; 30.03.2016