如何编辑 predict.lm() 使用的模型矩阵中的交互？
How to edit interactions in model matrix used by predict.lm()?
我想编辑 predict.lm() 在 R 中使用的模型矩阵来预测主效应而不是交互作用（但使用包含交互作用的完整模型的系数和方差）。
我试过：
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")
... 但是这个 returns 错误 '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")
使用@jay.sf 的回答，我还设法创建了一个版本，如果模型中存在多个级别的因素，我也可以使用该版本：
##full model (using block as a multi-level factor):
data(npk)
mod1 <- lm(yield ~ N*block, data=npk, x=T)
## get model formula and use it to generate the model matrix:
predgrid <- data.frame(model.matrix(mod1, data=npk))
## make a new dataframe using the model matrix and the response,
## and run the model using all columns in the new dataframe as terms:
npk2 <- as.data.frame(cbind(npk$yield, predgrid[, -1]))
colnames(npk2)[1] <- "yield"
mod2 <- lm(yield~., data=npk2)
## extract the model matrix dataframe again, to modify for predictions:
newmat <- predgrid[, -1]
colnames(newmat)
newmat[, 7:11] <- 0
## predict on modified matrix dataframe:
pred <- predict(mod2, newdata=newmat, type="response", interval="confidence")
head(pred) ##
# fit lwr upr
#1 48.15 41.18475 55.11525
#2 59.90 52.93475 66.86525
#3 48.15 41.18475 55.11525
#4 59.90 52.93475 66.86525
#5 67.50 55.43584 79.56416
#6 67.50 55.43584 79.56416
我想编辑 predict.lm() 在 R 中使用的模型矩阵来预测主效应而不是交互作用（但使用包含交互作用的完整模型的系数和方差）。
我试过：
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")
... 但是这个 returns 错误 '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")
使用@jay.sf 的回答，我还设法创建了一个版本，如果模型中存在多个级别的因素，我也可以使用该版本：
##full model (using block as a multi-level factor):
data(npk)
mod1 <- lm(yield ~ N*block, data=npk, x=T)
## get model formula and use it to generate the model matrix:
predgrid <- data.frame(model.matrix(mod1, data=npk))
## make a new dataframe using the model matrix and the response,
## and run the model using all columns in the new dataframe as terms:
npk2 <- as.data.frame(cbind(npk$yield, predgrid[, -1]))
colnames(npk2)[1] <- "yield"
mod2 <- lm(yield~., data=npk2)
## extract the model matrix dataframe again, to modify for predictions:
newmat <- predgrid[, -1]
colnames(newmat)
newmat[, 7:11] <- 0
## predict on modified matrix dataframe:
pred <- predict(mod2, newdata=newmat, type="response", interval="confidence")
head(pred) ##
# fit lwr upr
#1 48.15 41.18475 55.11525
#2 59.90 52.93475 66.86525
#3 48.15 41.18475 55.11525
#4 59.90 52.93475 66.86525
#5 67.50 55.43584 79.56416
#6 67.50 55.43584 79.56416