如何编辑 predict.lm() 使用的模型矩阵中的交互作用？

玉米旗布

我想编辑 R 中 predict.lm() 使用的模型矩阵来预测主效应而不是交互作用（但使用包含交互作用的完整模型的系数和方差）。

我试过了：

data(npk) #example data
mod <- lm(yield ~ N*P*K, data=npk, x=T) #run model
newmat <- mod$x # acquire model matrix
newmat[, c(5:8)] <- 0 #set interaction terms to 0
#try to predict on the new matrix..
predict(mod, as.data.frame(newmat), type="response", interval="confidence")

...但这会返回错误，'data' must be a data.frame, not a matrix or an array因为 predict.lm() 不接受模型矩阵。

如何使用示例代码中给出的模型矩阵进行预测？

（或者是否有更好的方法来预测主效应而不是相互作用，使用完整模型yield ~ N*P*K?）

周杰伦

我们可以手工计算交互作用；通过首先创建术语trms，然后以一种eval(parse())方法评估它们，可以轻松完成。

## create interaction terms 
iv <- c('N', 'P', 'K')  ## indp. vars
trms <- unlist(sapply(2:3, function(m) combn(iv, m, FUN=paste, collapse='x')))
## evaluate them to a matrix
Ia <- with(npk1, sapply(trms, function(x) eval(parse(text=gsub('x', '*', x)))))

然后只需 cbind 并在中使用它lm()，比较：

## cbind
npk2 <- cbind(npk1, Ia)
## following yield the same:
(mod1 <- lm(yield ~ .^3, data=npk1))
(mod2 <- lm(yield ~ ., data=npk2, x=TRUE))

那么你可以按照你的方法：

newmat <- mod2$x  ## acquire model matrix
newmat[, c(5:8)] <- 0  ## set interaction terms to 0
predict(mod2, newdata=as.data.frame(newmat))  ## newdata w/ Ia to zero
# 1        2        3        4        5        6        7        8        9       10 
# 54.90000 66.66667 51.43333 64.33333 63.76667 67.23333 52.00000 54.33333 54.33333 67.23333 
# 11       12       13       14       15       16       17       18       19       20 
# 63.76667 52.00000 63.76667 67.23333 52.00000 54.33333 66.66667 51.43333 64.33333 54.90000 
# 21       22       23       24 
# 64.33333 66.66667 54.90000 51.43333

然而：

predict(mod1)  ## old model
# 1        2        3        4        5        6        7        8        9       10 
# 50.50000 57.93333 51.43333 54.66667 63.76667 54.36667 52.00000 54.33333 54.33333 54.36667 
# 11       12       13       14       15       16       17       18       19       20 
# 63.76667 52.00000 63.76667 54.36667 52.00000 54.33333 57.93333 51.43333 54.66667 50.50000 
# 21       22       23       24 
# 54.66667 57.93333 50.50000 51.43333

数据：

npk1 <- structure(list(N = c(0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 
0, 0, 1, 0, 1, 0, 1, 1, 0, 0), P = c(1, 1, 0, 0, 0, 1, 0, 1, 
1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0), K = c(1, 0, 
0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 
0), yield = c(49.5, 62.8, 46.8, 57, 59.8, 58.5, 55.5, 56, 62.8, 
55.8, 69.5, 55, 62, 48.8, 45.5, 44.2, 52, 51.5, 49.8, 48.8, 57.2, 
59, 53.2, 56)), row.names = c(NA, 24L), class = "data.frame")

本文收集自互联网，转载请注明来源。

如有侵权，请联系 [email protected] 删除。