I am attempting to fit a model with a large number of predictors, such that it would be tedious to enumerate them in a model formula. This is straightforward to do with lm():

indicatorMatrix <- data.frame(matrix(rbinom(26000, 1, 1/3), ncol = 26))

colnames(indicatorMatrix) <- LETTERS

someDV <- rnorm(nrow(indicatorMatrix))

head(indicatorMatrix)

# One method, enumerating variables by name:

olsModel1 <- lm(someDV ~ A + B + C + D,  # ...etc.

               data = indicatorMatrix)

# Preferred method, including the matrix of predictors:

olsModel2 <- lm(someDV ~ as.matrix(indicatorMatrix))

summary(olsModel2)

Since I have a very large number of predictors (more than the 26 in this invented example), I don't want to list them individually as in the first example (someDV ~ A + B + C + D...), and I can avoid this by just including the predictors as.matrix.

However, I want to fit a mixed effects model, like this:

library(lme4)

meModel1 <- lmer(someDV ~ (1 | A) + (1 | B) + (1 | C),  # ...etc.

                 data = indicatorMatrix)

summary(meModel1)

Except that I want to include a large number of random effects terms. Rather than having to type (1 | A) ... (1 | ZZZ), I would like to include each predictor in a manner analogous to the matrix approach used for olsModel2 above. The following, obviously, does not work:

meModel2 <- lmer(someDV ~ (1 | as.matrix(indicatorMatrix)))

Do you have any suggestions for how I can best replicate the matrix-predictor approach for random effects with lmer()? I am very willing to consider "pragmatic" solutions (i.e. hacks), so long as they are "programmatic," and don't require me to copy & paste, etc. etc.

Thanks in advance for your time.