R中的分层多项式回归答案

【问题标题】：Hierarchical polynomial regression in RR中的分层多项式回归
【发布时间】：2015-04-23 13:32:24
【问题描述】：

我一直在尝试在 R 中拟合 sequential 多项式回归模型，但遇到了以下问题：虽然poly(x) 提供了一种快速方法，但该函数不遵守分层原则，大致指出在移动到高阶之前，所有低阶项都应该包含在模型中。

这里有一个解决方案，可能是自己手动选择进入模型的顺序，就像我对下面的玩具数据集所做的那样

pred<-matrix(c(rnorm(30),rnorm(30)),ncol=2)
y<-rnorm(30)

polys<-poly(pred,degree=4,raw=T)
z<-matrix(c(
#order 2
polys[,2],polys[,6],polys[,9],
#order 3
polys[,3],polys[,7],polys[,10],polys[,12],
#order 4
polys[,4],polys[,8],polys[,11],polys[,13],polys[,14]),
ncol=12)

polyreg3<-function(x){
BICm<-rep(0,dim(x)[2])
for(i in 1:dim(x)[2]){
model<-lm(y~pred[,1]+pred[,2]+x[,1:i]) #include one additional term each time
BICm[i]<-BIC(model)
}
list(BICm=BICm)
}

polyreg3(z)
which.min(polyreg3(z)$BICm)

但这对于较大次数的多项式在很大程度上是不切实际的。我当时想知道，有没有办法解决这个问题，最好是通过修改我的代码？

【问题讨论】：

for 循环最好在R 中避免。删除你的循环将是一回事。有很多关于如何在 SO 上执行此操作的示例（例如 here's a more generic example 或 one where somebody is applying a lm to a data.frame。此外，您可能希望使用 profr package 分析您的代码以找到您的瓶颈。
@RichardErickson 感谢您的建议，尽管它们不是我目前最紧迫的问题。

标签： r statistics

【解决方案1】：

如果我理解正确，您不仅需要原始自变量，还需要可以在给定学位的情况下创建的所有变量组合。

数据除以三个因变量、原始自变量和由model.frame() 创建的额外变量，给定一个度数（为简单起见，这里是 2）。

然后所有额外变量的组合都由combn() 和Map() 获得，因为选择列的方式是可变的（1 到 # 列）。

要拟合的数据集由cbind()创建，它们是自变量（ind）和原始自变量（original）和额外的变量组合（额外的）。

最后lm() 是合适的，得到BIC() 值。

如果目标是更高级别，则需要多次试验。例如，如果度数为 3，则应同时应用 2 度和 3 度。

set.seed(1237)
# independent variable
des <- data.frame(y = rnorm(30))
# dependent variables
pred<-matrix(c(rnorm(30), rnorm(30)), ncol=2)
# model frame given degree, 4095 combinations when degree = 4, set degree = 2 for simplicity
polys <- as.data.frame(poly(pred, degree = 2, raw = T))
# original independent variables
original <- polys[,c(names(polys)[names(polys) == "1.0" | names(polys) == "0.1"])]
# extra variables made by model.frame()
extra <- polys[,c(names(polys)[names(polys) != "1.0" & names(polys) != "0.1"])]
# all combinations of extra variables
# Map() for variable q in nCq, do.call() to make list neat
com <- do.call(c, Map(combn, ncol(extra), 1:ncol(extra), simplify = FALSE))
com
[[1]]
[1] 1

[[2]]
[1] 2

[[3]]
[1] 3

[[4]]
[1] 1 2

[[5]]
[1] 1 3

[[6]]
[1] 2 3

[[7]]
[1] 1 2 3

# data combined, followed by fitting lm()
bic <- lapply(com, function(x) {
  data <- cbind(des, original, extra[, x, drop = FALSE])
  BIC(lm(y ~ ., data))
})

do.call(c, bic)
[1] 100.3057 104.6485 104.8768 103.6572 103.4162 108.0270 106.7262

【讨论】：