Формула с динамическим числом переменных

Предположим, есть некоторый data.frame foo_data_frame , и нужно найти регрессию целевого столбца Y по некоторым другим столбцам. Для этого обычно используются некие формулы и модели. Например:

linear_model <- lm(Y ~ FACTOR_NAME_1 + FACTOR_NAME_2, foo_data_frame)

Это хорошо работает, если формула кодируется статически. Если требуется выполнить рутирование по нескольким моделям с постоянным числом зависимых переменных (скажем, 2), это можно рассматривать следующим образом:

for (i in seq_len(factor_number)) {
  for (j in seq(i + 1, factor_number)) {
    linear_model <- lm(Y ~ F1 + F2, list(Y=foo_data_frame$Y,
                                         F1=foo_data_frame[[i]],
                                         F2=foo_data_frame[[j]]))
    # linear_model further analyzing...
  }
}

Мой вопрос заключается в том, как сделать то же самое, когда количество переменных изменяется динамически во время программы. работает?

for (number_of_factors in seq_len(5)) {
   # Then root over subsets with #number_of_factors cardinality.
   for (factors_subset in all_subsets_with_fixed_cardinality) {
     # Here I want to fit model with factors from factors_subset.
     linear_model <- lm(Does R provide smth to write here?)
   }
}