多项式和交互项的岭回归与 glmnet答案

【问题标题】：Ridge Regression with glmnet for polynomial and interactions terms多项式和交互项的岭回归与 glmnet
【发布时间】：2018-12-03 20:46:18
【问题描述】：

我有一个包含 9000 个观察值和 13 个解释变量的数据集。

其中一些是分类变量，因此我将它们转换为虚拟变量并始终将一个类别设置为 NULL，因为它们是基本类别。现在我已经有了 53 个解释变量。我想做一个岭回归来获得样本外预测的最佳模型。为此，我想使用 glmnet 包。在 13 个解释变量中，我想创建 2-10 阶的多项式，并构建正常变量和所有多项式变量的所有可能交互项。我还想要 2-10 阶交互项的多项式。

我的问题是glmnet 包只使用矩阵或数据框作为参数，所以我不能使用公式。如果我尝试生成包含所有这些变量的数据框，我的数据框中会出现很多列，以至于我的 R 会关闭。

我能做些什么来解决这个问题？

【问题讨论】：

标签： r dataframe machine-learning regression

【解决方案1】：

您可以使用model.matrix 来指定它。这是一个使用iris数据的例子

首先，我将创建一个包含一些级别的虚拟因子列：

df <- iris
df$factor <- as.factor(sample(1:2, nrow(iris), replace = TRUE))
head(df)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species factor
1          5.1         3.5          1.4         0.2  setosa      1
2          4.9         3.0          1.4         0.2  setosa      1
3          4.7         3.2          1.3         0.2  setosa      1
4          4.6         3.1          1.5         0.2  setosa      2
5          5.0         3.6          1.4         0.2  setosa      2
6          5.4         3.9          1.7         0.4  setosa      2

现在制作模型矩阵

x <- model.matrix(Species ~ poly(Sepal.Length, 10)*factor-1, #-1 means no intercept
              data = df)
ncol(x) #output is 22, do head(x) to see what the columns are

基本上你有 10 列用于 Sepal.Length 的 poly 与第一级交互 + 10 列用于 Sepal.Length 的 poly 与第二级交互 + 两个对应于因子变量的热列。

现在使用 cv 找到最佳 lambda：

model <- cv.glmnet(y = iris$Species,
                   x = x,
                   alpha = 0,
                   family = "multinomial",
                   lambda.min.ratio = 1e-6) #changed it from the default since it looked the optimum is lower then the min lambda tried

plot(model)

【讨论】：

我已经用虚拟包创建了虚拟变量。我还创建了多项式并 cbind 两个数据集。我现在有一个包含 70 个 collums 的数据框，并且想要构建所有可能的交互项。模型矩阵不起作用，因为结果会大于 1 gb
你试过sparse.model.matrix吗？

【解决方案2】：

这不是一个完整的答案 - 更多的是解释为什么这如此困难 - 但肯定太长了，无法发表评论。

首先，lasso/ridge/elastic net 方法基本上建立在模型矩阵上的操作（即每个观察一行和一个 predictor 变量的矩阵，即从输入变量），因此在某些时候确实无法构建模型矩阵（尽管您可能可以分段进行，见下文）。

稀疏模型矩阵模型将有助于构建主要涉及因子（被转换为指标或虚拟变量）的模型矩阵，但我认为它不会对您有太大帮助。

这个math overflow question 解释了k 变量的顺序多项式n 需要choose(n+k,k) 派生变量（在R 术语中：这是给出k 来自@987654328 的可能样本数量的二项式系数@ 对象）。构造一个函数，报告模型矩阵的列数、元素总数（如果有r 行）以及模型矩阵的总大小（以 Mb 为单位）：

calc_size <- function(deg,nvars,r=9000) {
  cc <- choose(deg+nvars,nvars)
  return(c(cc,cc*r,cc*r*8/2^20))
}
calc_size(10,13,r=9000)

告诉您，您需要 1.14 百万列、102 十亿个总条目和 77 Gb 的空间来存储此问题的模型矩阵。这无需担心您需要多少空间来扩展分类变量。如果您真的想将 53 个全数字列处理到 10 阶，则需要 1.27*10^11 列和 PB 的存储空间。（稀疏模型矩阵较少，但组合稀疏（虚拟变量）和非稀疏列可能会很棘手......）

如果你真的想这样做，你可以使用biglasso 包。 vignette 给出了一个从文件支持的 31 G 数据集拟合 n = 2898, p = 1,339,511 的数据集的示例。这比你预期的要小，但至少它是相同的数量级（在 4 个核心上拟合大约需要 51 分钟......）如果我这样做，我会首先构建我的模型矩阵（可能是小块，例如一次 500 或 1000 行）并将块存储/连接到磁盘上的数据文件中，然后使用 biglasso 来拟合模型。

不过，根据您可用的硬件和技术专长水平，您可能不得不缩减您的目标（C(10,53) 列上的蛮力不太可能奏效）。

【讨论】：