r - 当保留高阶参数时，如何删除模型中的低阶参数？

Question

问题:只要高阶参数(即交互)保留在模型中，我就无法删除模型中的低阶参数(例如，主要效果参数)。即使这样做，模型也会被重构，并且新模型不会嵌套在更高的模型中。
请参见以下示例(由于我来自ANOVA，因此我使用contr.sum):

d <- data.frame(A = rep(c("a1", "a2"), each = 50), B = c("b1", "b2"), value = rnorm(100))
options(contrasts=c('contr.sum','contr.poly'))
m1 <- lm(value ~ A * B, data = d)
m1

## Call:
## lm(formula = value ~ A * B, data = d)
## 
## Coefficients:
## (Intercept)           A1           B1        A1:B1  
##   -0.005645    -0.160379    -0.163848     0.035523  

m2 <- update(m1, .~. - A)
m2

## Call:
## lm(formula = value ~ B + A:B, data = d)

## Coefficients:
## (Intercept)           B1       Bb1:A1       Bb2:A1  
##   -0.005645    -0.163848    -0.124855    -0.195902  

可以看出，尽管我删除了一个参数( A)，但新模型( m2)被重构，并且 并未将嵌套在更大的模型( m1)中。如果我将自己的因子转换为数值对比变量，则可以得到理想的结果，但是如何使用R的因子能力得到它呢？

问题:我如何删除R中的低阶因数并获得真正错过此参数且未重构的模型(即，较小模型中的参数数量必须更少)？

但是为什么呢？ 我想使用 lmer包中的 KRmodcomp函数为 pbkrtest模型获取类似于p值的“类型3”。所以这个例子实际上只是一个例子。

为什么不进行交叉验证？ 我觉得这实际上更像是一个R，而不是一个统计问题(即，我知道您永远不应该适合具有交互作用但没有主要影响的模型，但我仍然想这样做)。

Answer 1

这是一个答案。我所知道的没有办法直接通过公式来公式化该模型...

如上构造数据:

d <- data.frame(A = rep(c("a1", "a2"), each = 50),
                B = c("b1", "b2"), value = rnorm(100))
options(contrasts=c('contr.sum','contr.poly'))

确认原始发现，仅从公式中减去因子不起作用:

m1 <- lm(value ~ A * B, data = d)
coef(m1)
## (Intercept)          A1          B1       A1:B1 
## -0.23766309  0.04651298 -0.13019317 -0.06421580 

m2 <- update(m1, .~. - A)
coef(m2)
## (Intercept)          B1      Bb1:A1      Bb2:A1 
## -0.23766309 -0.13019317 -0.01770282  0.11072877 

制定新的模型矩阵:

X0 <- model.matrix(m1)
## drop Intercept column *and* A from model matrix
X1 <- X0[,!colnames(X0) %in% "A1"]

lm.fit允许直接指定模型矩阵:

m3 <- lm.fit(x=X1,y=d$value)
coef(m3)
## (Intercept)          B1       A1:B1 
## -0.2376631  -0.1301932  -0.0642158 

此方法仅适用于允许显式指定模型矩阵的一些特殊情况(例如 lm.fit， glm.fit)。

更普遍:

## need to drop intercept column (or use -1 in the formula)
X1 <- X1[,!colnames(X1) %in% "(Intercept)"]
## : will confuse things -- substitute something inert
colnames(X1) <- gsub(":","_int_",colnames(X1))
newf <- reformulate(colnames(X1),response="value")
m4 <- lm(newf,data=data.frame(value=d$value,X1))
coef(m4)
## (Intercept)          B1   A1_int_B1 
##  -0.2376631  -0.1301932  -0.0642158 

这种方法的缺点是，它不能将多个输入变量识别为来自相同的预测变量(即，多于2个水平因子的多个因子水平)。