我想自动确定 lm
中的哪些系数属于一个因子。所以假设我有以下模型:
d <- data.frame(a = gl(4, 2, 16), b = gl(2, 1, 16),
x = runif(16), y = runif(16), Y = runif(16))
l1 <- lm(Y ~ a + b + x + y, data = d)
l2 <- lm(Y ~ x + y, data = d)
那么第一个模型的系数名称如下:
names(coef(l1))
# [1] "(Intercept)" "a2" "a3" "a4" "b2"
# [6] "x" "y"
现在理想情况下,我想要一个函数告诉我 a2、a3、a4
和 b2
是虚拟编码因子的系数。
对于不包含任何因子的模型(如 l2
),输出应为 NULL
。
我查看了 str(l1)
,我发现有一个插槽 xlevels
(如果模型中存在因子)。我可以使用 names(l1$xlevels)
获取模型中所有因素的列表,然后在系数名称上使用 grep
:
names(coef(l1))[unlist(sapply(names(l1$xlevels), function(.) grep(., names(coef(l1)))))]
# [1] "a2" "a3" "a4" "b2"
但在我看来,这似乎是一个非常肮脏的解决方法,一旦我的模型中有相似的名称,它就不会起作用:
d$a4 <- runif(16)
l3 <- update(l1, . ~ . + a4, data = d)
names(coef(l3))[unlist(sapply(names(l3$xlevels), function(.) grep(., names(coef(l3)))))]
# [1] "a2" "a3" "a4" "a4" "b2"
此外,更改默认对比度会更改模型中虚拟系数的名称,因此即使是处理系数名称的最精细的策略也可能行不通。
长话短说:我如何获得属于一个因子的所有系数的列表?
最佳答案
这里有一些方法:
1) 这假定 model.matrix 中仅包含零和一的任何列都属于一个因子(截距除外)。它适用于 l1
、l2
和 l3
,非常短,不依赖于名称(拦截除外)并且不需要摆弄与 lm
对象组件。它适用于主效应和交互作用,因为如果主效应为 0/1,则交互作用也将如此。注释中的 l4
是 0/1 假设不成立的示例。
m <- model.matrix(l1)
all01 <- apply(m == 0 | m == 1, 2, all)
setdiff(names(all01[all01]), "(Intercept)")
## [1] "a2" "a3" "a4" "b2"
2) 这不使用名称(拦截除外)并且适用于 l1
、l2
和 l3
(和评论中的 l4
)。它不对模型矩阵做出任何假设,但仅适用于仅具有主要效果的模型。 (无拦截情况未经测试。)
cls <- attr(terms(l1), "dataClass")
intercept <- if ("(Intercept)" %in% names(coef(l1))) "" else "+ 0"
fn <- function(nm) names(coef(update(l1, paste(". ~", nm, intercept))))
setdiff(unlist(lapply(names(cls)[cls == "factor"], fn)), "(Intercept)")
关于r - 获取属于一个因子的所有系数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/36143679/