Ридж Регрессия с glmnet для полиномиальных и взаимодействующих членов

Question

У меня есть набор данных с 9000 наблюдениями и 13 пояснительными переменными.

Некоторые из них являются категориальными переменными, поэтому я преобразовал их в пустышки и всегда устанавливал одну категорию NULL, потому что это базовая категория.Теперь у меня уже есть 53 объясняющие переменные.Я хочу провести регрессию гребня, чтобы получить лучшую модель для прогнозов вне выборки.Для этого я хочу использовать пакет glmnet.Из 13 объясняющих переменных я хочу создать полиномы до степеней 2-10 и построить все возможные члены взаимодействия нормальных переменных, а также всех полиномиальных переменных.Мне также нужны полиномы для членов взаимодействия до степеней 2-10.

Моя проблема в том, что пакет glmnet использует только матрицы или фреймы данных в качестве аргументов, поэтому я не могу использовать формулы.Если я пытаюсь создать фрейм данных со всеми этими переменными, в моем фрейме данных появляется так много столбцов, что мой R выключается.

Что я могу сделать, чтобы решить эту проблему?

Answer 1

Это не совсем ответ - скорее объяснение того, почему это так сложно - но определенно слишком долго для комментария.

Прежде всего, методы лассо / риджа / эластичной сети основаны на операциях над матрицей модели (т. Е. Матрица с одной строкой на наблюдение и одной на предиктор переменная, т.е. переменные, полученные из входных данных переменные), так что в какой-то момент нет никакого способа построить матрицу модели (хотя вы могли бы сделать это кусочно, см. ниже).

Разреженные модельные матричные модели помогут в построении матричных моделей, которые в основном включают факторы (которые преобразуются в индикаторные или фиктивные переменные), но я думаю, что это не слишком вам поможет.

Этот вопрос о математическом переполнении объясняет, что для полинома n из k переменных требуется choose(n+k,k) производных переменных (в терминах R: это биномиальный коэффициент, дающий количество возможных выборок размера k из n+k объектов). Построение функции, которая сообщает количество столбцов матрицы модели, общее количество элементов при наличии r строк и общий размер (в Мб) матрицы модели:

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 миллион столбцов, 10,2 миллиардов записей и 77 ГБ места для хранения матрицы модели для этой проблемы. Это не беспокоясь о том, сколько еще места нужно для расширения категориальных переменных. Если вы действительно хотите обработать 53 полностью числовых столбца до 10-го порядка, вам потребуется 1,27 * 10 ^ 11 столбцов и петабайты дискового пространства. (Меньше для разреженных матриц моделей, но комбинирование разреженных (фиктивных переменных) и не разреженных столбцов может быть сложным ...)

Если вы действительно хотите это сделать, вы можете использовать пакет biglasso . Виньетка дает пример подгонки набора данных с n = 2898, p = 1,339,511 из набора данных 31 G на основе файла. Это меньше, чем планировалось, но, по крайней мере, оно того же порядка (подгонка занимает около 51 минуты на 4 ядрах ...). Если бы я это делал, я бы сначала построил матрицу модели (возможно, небольшими кусками скажем, 500 или 1000 строк за раз) и сохраните / объедините фрагменты в файле данных на диске, затем используйте biglasso, чтобы соответствовать модели.

Однако, в зависимости от имеющегося у вас оборудования и уровня технических знаний, вам, возможно, придется свернуть свои цели (грубая сила в столбцах C (10,53) вряд ли сработает).

Answer 2

Вы можете использовать model.matrix, чтобы указать это. Вот пример использования iris data

Сначала я создам столбец фиктивного фактора с несколькими уровнями:

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. Длина взаимодействует с первым уровнем + 10 столбцов поли для Sepal. Длина взаимодействует со вторым уровнем + два горячих столбца, соответствующих факторной переменной.

и теперь используйте cv, чтобы найти лучшую лямбду:

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)