为什么我得到 NA 系数以及 `lm` 如何降低交互的参考水平
Why do I get NA coefficients and how does `lm` drop reference level for interaction
我想了解 R 如何确定线性模型中交互的参考组。考虑以下因素：
df <- structure(list(id = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L,
5L, 5L, 5L, 5L, 5L, 5L), .Label = c("1", "2", "3", "4", "5"), class = "factor"),
year = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L), .Label = c("1", "2"), class = "factor"),
treatment = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L), .Label = c("0", "1"), class = "factor"),
y = c(1.4068142116718, 2.67041187927052, 2.69166439169131,
3.56550324537293, 1.60021286173782, 4.26629963353237, 3.85741108250572,
5.15740731957689, 4.15629768365669, 6.14875441068499, 3.31277276551286,
3.47223277168367, 3.74152201649338, 4.02734382610191, 4.49388620764795,
5.6432833241724, 4.76639399631094, 4.16885857079297, 4.96830394378801,
5.6286092105837, 6.60521404151111, 5.51371821706176, 3.97244221149279,
5.68793413111161, 4.90457233598412, 6.02826151378941, 4.92468415350312,
8.23718422822134, 5.87695836962708, 7.47264895892575)), .Names = c("id",
"year", "treatment", "y"), row.names = c(NA, -30L), class = "data.frame")
lm(y ~ -1 + id + year + year:treatment, df)
#Coefficients:
# id1 id2 id3 id4
# 2.6585 3.9933 4.1161 5.3544
# id5 year2 year1:treatment1 year2:treatment1
# 6.1991 0.7149 -0.6317 NA
R 尝试估计完整的交互集，而不是始终忽略参考组。结果，我得到了 NA 的结果。
此外，R 与它删除的组不一致。我想估计一个在主效应和交互作用中具有相同省略组 (year1) 的模型。如何强制 R 从前面的模型中省略 year1 和 year1:treatment1？
我知道这个问题有几种解决方法（例如，手动创建所有变量并将它们写在模型的公式中）。但是我估计的实际模型是这个问题的复杂得多的版本，这样的解决方法会很麻烦。
R tries to estimate the full set of interactions instead of consistently omitting a reference group. As a result, I am getting NAs in the results.
这两者之间没有因果关系。你得到 NA 纯粹是因为你的变量有嵌套。
R is inconsistent with which groups it drops. I would like to estimate a model with the same omitted group (year1) in the main effect and interactions. How to force R to omit year1 and year1:treatment1 from the preceding model?
没有不一致，但是model.matrix有自己的规律。你得到看似“不一致”的对比，因为你没有主效应 treatment 但只有交互作用 treatment:year.
下面我会先解释NA个系数，然后是交互的对比，最后给出一些建议。
NA 系数
默认对对比因子进行对比处理，默认contr.treatement，第一个因子水平作为参考水平。查看您的数据：
levels(df$id)
# [1] "1" "2" "3" "4" "5"
levels(df$year)
# [1] "1" "2"
levels(df$treatment)
# [1] "0" "1"
现在来看一个简单的线性模型：
lm(y ~ id + year + treatment, df)
#Coefficients:
#(Intercept) id2 id3 id4 id5 year2
# 2.153 1.651 1.773 2.696 3.541 1.094
# treatment1
# NA
可以看到id1、year1、treatment0没有，只是作为参考。
即使没有互动，你也已经有了一个NA系数。这意味着 treatment 嵌套在 span{id, year} 中。当您进一步包含 treatment:year 时，这样的嵌套仍然存在。
其实进一步测试发现treatment嵌套在id中：
lm(y ~ id + treatment, df)
# Coefficients:
#(Intercept) id2 id3 id4 id5 treatment1
# 2.700 1.651 1.773 2.696 3.541 NA
总而言之，变量 treatment 对于您的建模目的来说是完全多余的。如果包含 id，则无需包含 treatment 或 treatment:*，其中 * 可以是任何变量。
当回归模型中有很多因子变量时，很容易出现嵌套。注意，对比不一定会消除嵌套，因为对比只识别变量名，而不是潜在的数字特征。 看下面的例子如何作弊contr.treatment：
df$ID <- df$id
lm(y ~ id + ID, df)
#Coefficients:
#(Intercept) id2 id3 id4 id5 ID2
# 2.700 1.651 1.773 2.696 3.541 NA
# ID3 ID4 ID5
# NA NA NA
看，对比度按预期工作，但是 ID 嵌套在 id 中，所以我们有排名不足。
交互对比
我们首先通过删除 id 变量来消除 NA 施加的噪音。然后，具有 treatment 和 year 的回归模型将是满秩的，因此如果对比成功，则不应看到 NA。
相互作用或高阶效应的对比取决于低阶效应的存在。比较以下模型：
lm(y ~ year:treatment, df) ## no low-order effects at all
#Coefficients:
# (Intercept) treatment0:year1 treatment1:year1 treatment0:year2
# 5.4523 -1.3976 -1.3466 -0.6826
#treatment1:year2
# NA
lm(y ~ year + treatment + year:treatment, df) ## with all low-order effects
#Coefficients:
# (Intercept) treatment1 year2 treatment1:year2
# 4.05471 0.05094 0.71493 0.63170
在第一个模型中，没有做对比，因为没有主效应，只有二阶效应。这里的NA是因为没有对比。
现在考虑另一组例子，包括一些但不是所有的主要影响：
lm(y ~ year + year:treatment, df) ## main effect `treatment` is missing
#Coefficients:
# (Intercept) year2 year1:treatment1 year2:treatment1
# 4.05471 0.71493 0.05094 0.68264
lm(y ~ treatment + year:treatment, df) ## main effect `year` is missing
#Coefficients:
# (Intercept) treatment1 treatment0:year2 treatment1:year2
# 4.05471 0.05094 0.71493 1.34663
在这里，交互作用的对比将发生在主效应缺失的变量上。例如，在第一个模型中，主效应 treatment 缺失，因此交互作用下降 treatement0；而在第二个模型中，主效应 year 缺失，因此交互作用下降 year1.
一般准则
指定高阶效应时始终包括所有低阶效应。这不仅给出了易于理解的对比行为，而且还有一些其他吸引人的统计原因。您可以在 Cross Validated 上阅读 Including the interaction but not the main effects in a model。
另一个建议是始终包含拦截。在线性回归中，具有截距的模型会产生无偏估计，而残差的均值为 0。
我想了解 R 如何确定线性模型中交互的参考组。考虑以下因素：
df <- structure(list(id = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L,
5L, 5L, 5L, 5L, 5L, 5L), .Label = c("1", "2", "3", "4", "5"), class = "factor"),
year = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L), .Label = c("1", "2"), class = "factor"),
treatment = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L), .Label = c("0", "1"), class = "factor"),
y = c(1.4068142116718, 2.67041187927052, 2.69166439169131,
3.56550324537293, 1.60021286173782, 4.26629963353237, 3.85741108250572,
5.15740731957689, 4.15629768365669, 6.14875441068499, 3.31277276551286,
3.47223277168367, 3.74152201649338, 4.02734382610191, 4.49388620764795,
5.6432833241724, 4.76639399631094, 4.16885857079297, 4.96830394378801,
5.6286092105837, 6.60521404151111, 5.51371821706176, 3.97244221149279,
5.68793413111161, 4.90457233598412, 6.02826151378941, 4.92468415350312,
8.23718422822134, 5.87695836962708, 7.47264895892575)), .Names = c("id",
"year", "treatment", "y"), row.names = c(NA, -30L), class = "data.frame")
lm(y ~ -1 + id + year + year:treatment, df)
#Coefficients:
# id1 id2 id3 id4
# 2.6585 3.9933 4.1161 5.3544
# id5 year2 year1:treatment1 year2:treatment1
# 6.1991 0.7149 -0.6317 NA
R 尝试估计完整的交互集，而不是始终忽略参考组。结果，我得到了 NA 的结果。
此外，R 与它删除的组不一致。我想估计一个在主效应和交互作用中具有相同省略组 (year1) 的模型。如何强制 R 从前面的模型中省略 year1 和 year1:treatment1？
我知道这个问题有几种解决方法（例如，手动创建所有变量并将它们写在模型的公式中）。但是我估计的实际模型是这个问题的复杂得多的版本，这样的解决方法会很麻烦。
R tries to estimate the full set of interactions instead of consistently omitting a reference group. As a result, I am getting
NAs in the results.
这两者之间没有因果关系。你得到 NA 纯粹是因为你的变量有嵌套。
R is inconsistent with which groups it drops. I would like to estimate a model with the same omitted group (
year1) in the main effect and interactions. How to force R to omityear1andyear1:treatment1from the preceding model?
没有不一致，但是model.matrix有自己的规律。你得到看似“不一致”的对比，因为你没有主效应 treatment 但只有交互作用 treatment:year.
下面我会先解释NA个系数，然后是交互的对比，最后给出一些建议。
NA 系数
默认对对比因子进行对比处理，默认contr.treatement，第一个因子水平作为参考水平。查看您的数据：
levels(df$id)
# [1] "1" "2" "3" "4" "5"
levels(df$year)
# [1] "1" "2"
levels(df$treatment)
# [1] "0" "1"
现在来看一个简单的线性模型：
lm(y ~ id + year + treatment, df)
#Coefficients:
#(Intercept) id2 id3 id4 id5 year2
# 2.153 1.651 1.773 2.696 3.541 1.094
# treatment1
# NA
可以看到id1、year1、treatment0没有，只是作为参考。
即使没有互动，你也已经有了一个NA系数。这意味着 treatment 嵌套在 span{id, year} 中。当您进一步包含 treatment:year 时，这样的嵌套仍然存在。
其实进一步测试发现treatment嵌套在id中：
lm(y ~ id + treatment, df)
# Coefficients:
#(Intercept) id2 id3 id4 id5 treatment1
# 2.700 1.651 1.773 2.696 3.541 NA
总而言之，变量 treatment 对于您的建模目的来说是完全多余的。如果包含 id，则无需包含 treatment 或 treatment:*，其中 * 可以是任何变量。
当回归模型中有很多因子变量时，很容易出现嵌套。注意，对比不一定会消除嵌套，因为对比只识别变量名，而不是潜在的数字特征。 看下面的例子如何作弊contr.treatment：
df$ID <- df$id
lm(y ~ id + ID, df)
#Coefficients:
#(Intercept) id2 id3 id4 id5 ID2
# 2.700 1.651 1.773 2.696 3.541 NA
# ID3 ID4 ID5
# NA NA NA
看，对比度按预期工作，但是 ID 嵌套在 id 中，所以我们有排名不足。
交互对比
我们首先通过删除 id 变量来消除 NA 施加的噪音。然后，具有 treatment 和 year 的回归模型将是满秩的，因此如果对比成功，则不应看到 NA。
相互作用或高阶效应的对比取决于低阶效应的存在。比较以下模型：
lm(y ~ year:treatment, df) ## no low-order effects at all
#Coefficients:
# (Intercept) treatment0:year1 treatment1:year1 treatment0:year2
# 5.4523 -1.3976 -1.3466 -0.6826
#treatment1:year2
# NA
lm(y ~ year + treatment + year:treatment, df) ## with all low-order effects
#Coefficients:
# (Intercept) treatment1 year2 treatment1:year2
# 4.05471 0.05094 0.71493 0.63170
在第一个模型中，没有做对比，因为没有主效应，只有二阶效应。这里的NA是因为没有对比。
现在考虑另一组例子，包括一些但不是所有的主要影响：
lm(y ~ year + year:treatment, df) ## main effect `treatment` is missing
#Coefficients:
# (Intercept) year2 year1:treatment1 year2:treatment1
# 4.05471 0.71493 0.05094 0.68264
lm(y ~ treatment + year:treatment, df) ## main effect `year` is missing
#Coefficients:
# (Intercept) treatment1 treatment0:year2 treatment1:year2
# 4.05471 0.05094 0.71493 1.34663
在这里，交互作用的对比将发生在主效应缺失的变量上。例如，在第一个模型中，主效应 treatment 缺失，因此交互作用下降 treatement0；而在第二个模型中，主效应 year 缺失，因此交互作用下降 year1.
一般准则
指定高阶效应时始终包括所有低阶效应。这不仅给出了易于理解的对比行为，而且还有一些其他吸引人的统计原因。您可以在 Cross Validated 上阅读 Including the interaction but not the main effects in a model。
另一个建议是始终包含拦截。在线性回归中，具有截距的模型会产生无偏估计，而残差的均值为 0。