获取属于一个因子的所有系数

Question

我想自动判断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"

此外，更改默认对比度会更改我的模型中虚拟系数的名称，因此即使是处理系数名称的最详尽的策略也可能不起作用。

长话短说：我如何获得属于一个因子的所有系数的列表？

Answer 1

这里有一些方法：

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

Answer 2

经过评论区的讨论，我终于想出了这个解决方案。请注意，我将期望的结果略微更改为 return，不仅是分配给因子的系数，而且还区分了它们是属于因子主效应、因子-因子交互作用还是因子-变量交互作用。我将所有用例都包含在讨论之外，并且输出如预期的那样对系数进行了适当的表征。

getCoefficientType <- function(mod) {
   INTCPT <- "(Intercept)"
   te <- mod$terms
   hasIntercept <- attr(te, "intercept") == 1
   ## factor terms
   predictors <- attr(te, "dataClasses")
   factors <- names(predictors[predictors == "factor"])
   if (hasIntercept) {
      termLabels <- c(INTCPT, attr(te, "term.labels"))
   } else {
      termLabels <- attr(te, "term.labels")
   }
   ## - loop through all terms in the model
   ## - split interactions at ":" into atoms
   ## - check if any of the atoms occurs in the factor list
   types <- sapply(strsplit(termLabels, ":"), function(x) {
      ind <- x %in% factors
      if (length(x) == 1) {
         if (x == INTCPT) {
            "intercept"
         } else if (ind) {
            "factor.main"
         } else {
            "variable.main"
         }
      } else {
         if (all(ind)) {
            "factor.factor.interaction"
         } else if (!any(ind)) {
            "variable.variable.interaction"
         } else {
            "factor.variable.interaction"
         }
      }
   })
   setNames(rep(types, rle(mod$assign)$length), names(coef(mod)))
}

d <- data.frame(a = gl(4, 2, 16), b = gl(2, 1, 16), 
                x = runif(16), y = runif(16), Y = runif(16), a4 = runif(16))
l1 <- lm(Y ~ a + b + x + y, data = d)
l2 <- lm(Y ~ x + y, data = d)
l3 <- update(l1, . ~ . + a4, data = d)
l4 <- update(l3, contrasts = list(a = "contr.poly"))
l5 <- update(l2, . ~ . + a:x + x:y)
l6 <- update(l5, . ~ . - 1)
getCoefficientType(l1)
getCoefficientType(l2)
getCoefficientType(l3)
getCoefficientType(l4)
getCoefficientType(l5)
getCoefficientType(l6)

获取属于一个因子的所有系数

Get all Coefficients Belonging to a Factor

r

lm