确定哪个主效应属于哪个交互作用

Question

假设我们有以下模型：

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 关联的主效应，而不是因子 a 的虚拟编码因子 4。

因此，我展示的模型是 R 中的有效模型，但系数的名称不明确。那么有没有其他方法可以在交互系数和相应的主效应之间建立正确的映射？

Answer 1

您可以使用 lm 对象的 assign 组件：

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

这会将系数映射到扩展公式 a4 + a + x + a:x。

有关 assign 组件的文档，请参阅 help("model.matrix")。

编辑：

要扩展我的答案，您可以这样做：

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

确定哪个主效应属于哪个交互作用

Determine Which Main Effect Belongs to Which Interaction

interaction

r