如何在 R 中使用公式排除主效应但保留交互作用
How to use formula in R to exclude main effect but retain interaction
我不想要主效应，因为它与更精细的固定效应共线，所以有这些 NA.
很烦人
在这个例子中：
lm(y ~ x * z)
我想要 x（数字）和 z（因子）的交互作用，但不是 z 的主效应。
简介
?formula 的 R 文档说：
The ‘*’ operator denotes factor crossing: ‘a * b’ interpreted as ‘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 的变量 class。
- 如果其中none个是因数，或者只有一个是因数，则为真；
- 如果两者都是因数，则否。当存在交互作用（高阶效应）时，您永远不能放弃主效应（低阶效应）。
这种行为是由于一个名为 model.matrix.default 的神奇函数，它根据公式构造设计矩阵。数值变量仅包含在列中，但因子变量会自动编码为许多虚拟列。正是这种虚拟重新编码是一种魔法。人们普遍认为我们可以启用或禁用对比度来控制它，但事实并非如此。即使在 中，我们也会失去对对比的控制。问题是 model.matrix.default 在做虚拟编码时有自己的规则，它对你如何指定模型公式非常敏感。正是由于这个原因，当两个因素之间存在交互作用时，我们不能丢弃主效应。
数值和因子之间的相互作用
根据您的问题，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" 两个级别都是虚拟编码的。鉴于两个规范最终都具有三个系数，它们实际上是等价的（从数学上讲，两种情况下的设计矩阵具有相同的列space），您可以通过比较它们的拟合值来验证这一点：
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。它能够做出正确的决定！
两个因素之间的相互作用
我重申一下：没有办法在有交互的情况下降低主要效果。
我不会在这里提供额外的示例，因为我在 中有一个示例。请参阅那里的 "Contrasts 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 很烦人，因为它在虚拟编码中似乎没有一致的方式。他们认为 "consistent manner" 总是降低第一个因素水平。好吧，没问题，您可以通过手动进行虚拟编码来绕过 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
.
但是，我真的不推荐这种黑客攻击。当您将模型公式传递给 lm 等时，model.matrix.default 试图为您提供最合理的构造。此外，实际上我们想用拟合模型进行预测。如果您自己进行了虚拟编码，则将 newdata 提供给 predict.
时会遇到困难
这是一个很好的解释，但让我在选择重要预测变量的过程中再补充一件事。
让我们再次考虑以下模型：
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，众所周知，return 不是 class lm，而是一个list。所以问题是：你知道有什么方法可以在不丢失 lm class 的情况下消除主效应吗？
使用 lm.fit 但得到一个 lm class 对象：
我不是程序员，但 lm-function 的简单改编对我有用：我添加了 lm.fit 所需的两个参数（模型矩阵和响应变量）并替换了（在 lm-function) 我的模型矩阵 Xmat 和响应的 x 和 y（用于 lm 中的 lm.fit）：
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)
所有这一切当然可以更详细，但对我来说它有效 -> 它 returns 一个 lm-object 可以与 summary、resid、sim、plot ...
最佳
皮乌斯
我不想要主效应，因为它与更精细的固定效应共线，所以有这些 NA.
在这个例子中：
lm(y ~ x * z)
我想要 x（数字）和 z（因子）的交互作用，但不是 z 的主效应。
简介
?formula 的 R 文档说：
The ‘*’ operator denotes factor crossing: ‘a * b’ interpreted as ‘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 的变量 class。
- 如果其中none个是因数，或者只有一个是因数，则为真；
- 如果两者都是因数，则否。当存在交互作用（高阶效应）时，您永远不能放弃主效应（低阶效应）。
这种行为是由于一个名为 model.matrix.default 的神奇函数，它根据公式构造设计矩阵。数值变量仅包含在列中，但因子变量会自动编码为许多虚拟列。正是这种虚拟重新编码是一种魔法。人们普遍认为我们可以启用或禁用对比度来控制它，但事实并非如此。即使在 model.matrix.default 在做虚拟编码时有自己的规则，它对你如何指定模型公式非常敏感。正是由于这个原因，当两个因素之间存在交互作用时，我们不能丢弃主效应。
数值和因子之间的相互作用
根据您的问题，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" 两个级别都是虚拟编码的。鉴于两个规范最终都具有三个系数，它们实际上是等价的（从数学上讲，两种情况下的设计矩阵具有相同的列space），您可以通过比较它们的拟合值来验证这一点：
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。它能够做出正确的决定！
两个因素之间的相互作用
我重申一下：没有办法在有交互的情况下降低主要效果。
我不会在这里提供额外的示例，因为我在
~ 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 很烦人，因为它在虚拟编码中似乎没有一致的方式。他们认为 "consistent manner" 总是降低第一个因素水平。好吧，没问题，您可以通过手动进行虚拟编码来绕过 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
但是，我真的不推荐这种黑客攻击。当您将模型公式传递给 lm 等时，model.matrix.default 试图为您提供最合理的构造。此外，实际上我们想用拟合模型进行预测。如果您自己进行了虚拟编码，则将 newdata 提供给 predict.
这是一个很好的解释，但让我在选择重要预测变量的过程中再补充一件事。
让我们再次考虑以下模型：
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，众所周知，return 不是 class lm，而是一个list。所以问题是：你知道有什么方法可以在不丢失 lm class 的情况下消除主效应吗？
使用 lm.fit 但得到一个 lm class 对象： 我不是程序员，但 lm-function 的简单改编对我有用：我添加了 lm.fit 所需的两个参数（模型矩阵和响应变量）并替换了（在 lm-function) 我的模型矩阵 Xmat 和响应的 x 和 y（用于 lm 中的 lm.fit）：
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)
所有这一切当然可以更详细，但对我来说它有效 -> 它 returns 一个 lm-object 可以与 summary、resid、sim、plot ...
最佳 皮乌斯