如何在R中使用公式排除主效应但保留交互作用

Question

我不想要主效应，因为它与更精细的因子固定效应共线，所以拥有这些NA 很烦人。

在这个例子中：

lm(y ~ x * z)

我想要x（数字）和z（因子）的交互，但不是z的主要效果。

Answer 1

简介

?formula 的 R 文档说：

“*”运算符表示因子交叉：“a * b”解释为“a + b + a : b”

所以听起来删除主效应很简单，只需执行以下操作之一：

a + a:b  ## main effect on `b` is dropped
b + a:b  ## main effect on `a` is dropped
a:b      ## main effects on both `a` and `b` are dropped

哦，真的吗？不不不（太简单，太天真）。实际上它取决于a 和b 的变量类。

• 如果它们都不是因数，或者只有一个是因数，则为真；
• 如果它们都是因子，则不。当存在交互（高阶效应）时，您永远不能放弃主效应（低阶效应）。

这种行为是由名为model.matrix.default 的魔术函数引起的，该函数根据公式构造设计矩阵。数值变量只是按原样包含在列中，但因子变量会自动编码为许多虚拟列。正是这种虚拟的重新编码是一种魔法。通常认为我们可以启用或禁用对比度来控制它，但实际上并非如此。即使在this simplest example 中，我们也无法控制对比。问题是model.matrix.default在做dummy encoding时有自己的规则，对你如何指定模型公式非常敏感。正是由于这个原因，当两个因素之间存在交互时，我们不能丢弃主效应。

---

数字和因子之间的相互作用

根据您的问题，x 是数字，z 是一个因素。您可以通过

指定具有z主效应但具有交互作用的模型

y ~ x + x:z

由于x是数字，所以等价于do

y ~ x:z

这里唯一的区别是参数化（或者model.matrix.default 如何进行虚拟编码）。考虑一个小例子：

set.seed(0)
y <- rnorm(10)
x <- rnorm(10)
z <- gl(2, 5, labels = letters[1:2])

fit1 <- lm(y ~ x + x:z)
#Coefficients:
#(Intercept)            x         x:zb  
#     0.1989      -0.1627      -0.5456  

fit2 <- lm(y ~ x:z)
#Coefficients:
#(Intercept)         x:za         x:zb  
#     0.1989      -0.1627      -0.7082 

从系数的名称我们看到，在第一个规范中，z 是对比的，因此它的第一级“a”不是虚拟编码，而在第二个规范中，z 没有对比，两个级别“ a" 和 "b" 是虚拟编码的。鉴于这两个规范最终都有三个系数，它们实际上是等价的（从数学上讲，两种情况下的设计矩阵具有相同的列空间），您可以通过比较它们的拟合值来验证这一点：

all.equal(fit1$fitted, fit2$fitted)
# [1] TRUE

那么为什么z 在第一种情况下是相反的呢？因为否则我们有两个虚拟列x:z，这两列之和只是x，在公式中与现有模型项x 别名。事实上，在这种情况下，即使你要求不想要对比，model.matrix.default 也不会服从：

model.matrix.default(y ~ x + x:z,
      contrast.arg = list(z = contr.treatment(nlevels(z), contrasts = FALSE)))
#   (Intercept)          x       x:zb
#1            1  0.7635935  0.0000000
#2            1 -0.7990092  0.0000000
#3            1 -1.1476570  0.0000000
#4            1 -0.2894616  0.0000000
#5            1 -0.2992151  0.0000000
#6            1 -0.4115108 -0.4115108
#7            1  0.2522234  0.2522234
#8            1 -0.8919211 -0.8919211
#9            1  0.4356833  0.4356833
#10           1 -1.2375384 -1.2375384

那么为什么在第二种情况下z 没有对比？因为如果是这样，我们在构建交互时就会失去“a”级的效果。即使你需要对比，model.matrix.default 也会忽略你：

model.matrix.default(y ~ x:z,
      contrast.arg = list(z = contr.treatment(nlevels(z), contrasts = TRUE)))
#   (Intercept)       x:za       x:zb
#1            1  0.7635935  0.0000000
#2            1 -0.7990092  0.0000000
#3            1 -1.1476570  0.0000000
#4            1 -0.2894616  0.0000000
#5            1 -0.2992151  0.0000000
#6            1  0.0000000 -0.4115108
#7            1  0.0000000  0.2522234
#8            1  0.0000000 -0.8919211
#9            1  0.0000000  0.4356833
#10           1  0.0000000 -1.2375384

哦，太棒了model.matrix.default。它能够做出正确的决定！

---

两个因素之间的相互作用

让我重申一下：当存在交互时，没有办法放弃主效应。

我不会在这里提供额外的示例，因为我在Why do I get NA coefficients and how does lm drop reference level for interaction 中有一个示例。请参阅那里的“交互对比”部分。简而言之，以下所有规格都给出了相同的模型（它们具有相同的拟合值）：

~ year:treatment
~ year:treatment + 0
~ year + year:treatment
~ treatment + year:treatment
~ year + treatment + year:treatment
~ year * treatment

特别是，第一个规范导致NA 系数。

所以一旦~ 的RHS 包含year:treatment，你就永远不能要求model.matrix.default 放弃主效应。

People inexperienced with this behavior are to be surprised when producing ANOVA tables.

---

绕过model.matrix.default

有些人认为model.matrix.default 很烦人，因为它在虚拟编码中似乎没有一致的方式。在他们看来，“一致的方式”是始终降低第一个因素水平。好吧，没问题，您可以通过手动进行虚拟编码绕过model.matrix.default，并将生成的虚拟矩阵作为变量提供给lm等。

但是，您仍然需要model.matrix.default 的帮助才能轻松地对 a（是的，只有一个）因子变量进行虚拟编码。例如，对于我们之前示例中的变量z，其完整的虚拟编码如下，您可以保留其全部或部分列进行回归。

Z <- model.matrix.default(~ z + 0)  ## no contrasts (as there is no intercept)
#   za zb
#1   1  0
#2   1  0
#3   1  0
#4   1  0
#5   1  0
#6   0  1
#7   0  1
#8   0  1
#9   0  1
#10  0  1
#attr(,"assign")
#[1] 1 1
#attr(,"contrasts")
#attr(,"contrasts")$z
#[1] "contr.treatment"

回到我们的简单示例，如果我们不想在y ~ x + x:z 中对比z，我们可以这样做

Z2 <- Z[, 1:2]  ## use "[" to remove attributes of `Z`
lm(y ~ x + x:Z2)
#Coefficients:
#(Intercept)            x       x:Z2za       x:Z2zb  
#     0.1989      -0.7082       0.5456           NA

我们看到NA 并不奇怪（因为colSums(Z2) 与x 有别名）。如果我们想在y ~ x:z 中强制执行对比，我们可以执行以下任一操作：

Z1 <- Z[, 1]
lm(y ~ x:Z1)
#Coefficients:
#(Intercept)         x:Z1  
#    0.34728     -0.06571

Z1 <- Z[, 2]
lm(y ~ x:Z1)
#Coefficients:
#(Intercept)         x:Z1  
#     0.2318      -0.6860  

And the latter case is probably what contefranz is trying to do.

但是，我并不真正推荐这种黑客行为。当您将模型公式传递给lm 等时，model.matrix.default 试图为您提供最合理的构造。此外，实际上我们希望使用拟合模型进行预测。如果您自己完成了虚拟编码，则在将newdata 提供给predict 时会很困难。

Answer 2

这是一个很好的解释，但在选择重要预测变量的过程中，让我再补充一点。

让我们再次考虑以下模型：

fit1 <- lm(y ~ x + x:z)
#Coefficients:
#(Intercept)            x         x:zb  
#     0.1989      -0.1627      -0.5456

假设主效应x 在统计上不显着，而您想去掉它。至少对我来说，最直观的就是将第二个模型写成如上：

fit2 <- lm(y ~ x:z)
#Coefficients:
#(Intercept)         x:za         x:zb  
#     0.1989      -0.1627      -0.7082 

这最终将主效应恢复为与因子基线水平的交互作用。现在，为了真正不包括主效应，我能找到的唯一解决方案是利用lm.fit，众所周知，它不会返回类lm的对象，而是list。所以问题是：你知道有什么方法可以在不丢失lm 类的情况下摆脱主效应吗？

Answer 3

使用 lm.fit 但得到一个 lm 类对象： 我不是程序员，但对 lm 函数的简单改编对我有用：我添加了 lm.fit 所需的两个参数（模型矩阵和响应变量）并替换（在 lm 函数中）x 和y（用于 lm.fit 内的 lm）由我的模型矩阵 Xmat 和响应：

lm.2 <- function (formula, data, subset, weights, na.action, method = "qr", 
                  model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
                  contrasts = NULL, offset, ..., Xmat=NA, response=NA)
{
  ret.x <- x
  ret.y <- y
  cl <- match.call()
  mf <- match.call(expand.dots = FALSE)
  m <- match(c("formula", "data", "subset", "weights", "na.action", 
               "offset"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf$drop.unused.levels <- TRUE
  mf[[1L]] <- quote(stats::model.frame)
  mf <- eval(mf, parent.frame())
  if (method == "model.frame") 
    return(mf)
  else if (method != "qr") 
    warning(gettextf("method = '%s' is not supported. Using 'qr'", 
                     method), domain = NA)
  mt <- attr(mf, "terms")
  y <- model.response(mf, "numeric")
  w <- as.vector(model.weights(mf))
  if (!is.null(w) && !is.numeric(w)) 
    stop("'weights' must be a numeric vector")
  offset <- as.vector(model.offset(mf))
  if (!is.null(offset)) {
    if (length(offset) != NROW(y)) 
      stop(gettextf("number of offsets is %d, should equal %d (number of observations)", 
                    length(offset), NROW(y)), domain = NA)
  }
  if (is.empty.model(mt)) {
    x <- NULL
    z <- list(coefficients = if (is.matrix(y)) matrix(NA_real_, 
                                                      0, ncol(y)) else numeric(), residuals = y, fitted.values = 0 * 
                y, weights = w, rank = 0L, df.residual = if (!is.null(w)) sum(w != 
                                                                                0) else if (is.matrix(y)) nrow(y) else length(y))
    if (!is.null(offset)) {
      z$fitted.values <- offset
      z$residuals <- y - offset
    }
  }
  else {
    z <- if (is.null(w)) 
      lm.fit(Xmat, response, offset = offset, singular.ok = singular.ok, 
             ...)
    else lm.wfit(Xmat, response, w, offset = offset, singular.ok = singular.ok, 
                 ...)
  }
  class(z) <- c(if (is.matrix(y)) "mlm", "lm")
  z$na.action <- attr(mf, "na.action")
  z$offset <- offset
  z$contrasts <- attr(x, "contrasts")
  z$xlevels <- .getXlevels(mt, mf)
  z$call <- cl
  z$terms <- mt
  if (model) 
    z$model <- mf
  if (ret.x) 
    z$x <- x
  if (ret.y) 
    z$y <- y
  if (!qr) 
    z$qr <- NULL
  z
}

然后，照常创建模型矩阵，但删除不需要的列：

Xmat <- model.matrix(~x + x:z, data=mydata)
Xmat <- Xmat[,-(colnames(Xmat)=="x")]

mod <- lm.2(~x + x:z, data=mydata, Xmat=Xmat, response=mydata$y)

所有这些当然可以更详细，但对我来说它有效 -> 它返回一个 lm 对象，可以与 summary、resid、sim、plot 一起使用......

最好的 庇护