How to change combinations of random and fixed effects automatically for lmer formulae (in R)

Solution for How to change combinations of random and fixed effects automatically for lmer formulae (in R)
is Given Below:

I have previously been able to find the combinations of fixed effects,
See Here , and been able to incorporate random effects as needed, including nested random effects.

However, a colleague is convinced we need to check more models by investigating the variables as random effects. (Let’s not go into the validity of that here, as I am not that convinced, but I am prepared test it and check the resultant stats).
Currently I have this code:

first_col <- "taxon"
# The second column always needed, random effect 
second_col <- "(1|specimenN/A/B)"
other_cols <- c("C", "D", "E", "F")
secondLine <- paste(first_col, second_col, sep = " + ")

combos =, lapply(seq_along(other_cols), function(y) {
    arrangements::combinations(other_cols, y, layout = "l")
theBigList = sapply(combos, paste, collapse = " + ")
theBigList = paste(first_col, second_col, theBigList, sep = " + ")
# Once the combinations are done for the columns, we also need to add the ones that are at the start, but don't have combo's added

theBigList = c(second_col, secondLine, theBigList)
> theBigList
 [1] "(1|specimenN/A/B)"                         "taxon + (1|specimenN/A/B)"                 "taxon + (1|specimenN/A/B) + C"            
 [4] "taxon + (1|specimenN/A/B) + D"             "taxon + (1|specimenN/A/B) + E"             "taxon + (1|specimenN/A/B) + F"            
 [7] "taxon + (1|specimenN/A/B) + C + D"         "taxon + (1|specimenN/A/B) + C + E"         "taxon + (1|specimenN/A/B) + C + F"        
[10] "taxon + (1|specimenN/A/B) + D + E"         "taxon + (1|specimenN/A/B) + D + F"         "taxon + (1|specimenN/A/B) + E + F"        
[13] "taxon + (1|specimenN/A/B) + C + D + E"     "taxon + (1|specimenN/A/B) + C + D + F"     "taxon + (1|specimenN/A/B) + C + E + F"    
[16] "taxon + (1|specimenN/A/B) + D + E + F"     "taxon + (1|specimenN/A/B) + C + D + E + F"

theBigList is sent to a function that adds the dependent variable as it works through the list (but we can essentially ignore this aspect of it for these purposes, but it necessary to be like this to run through the many dependent variables loop) i.e.

modelName <- models[i]
model.formula <- paste(response, modelName, sep=' ~ ')
fit <- lmer(model.formula, data=df, REML=F)

I would like to be able to generate random effects for "A", "B", "C", "D", "E", "F" without them repeating as fixed effects i.e. A OR (1|A) etc.

Theoretically, the inputs could be something like this:

# The First column that is ALWAYS in the models tested. It won't be removed. Usually the main Fixed effect (i.e. variables of primary interest)
first_col <- "taxon"

# All possible random effects with nestings you want to try, that INCLUDE your PRIMARY random effect. One of these will always be in the model
nestings1 <- c("(1|specNum)", "(1|specNum/A/B)", "(1|specNum/A)")

#other possible "valid" nestings of random effects you want to try
nestings2 <- c("(1|A/B)", "(1|D/F)")

# All the Fixed Effects you would like to test. All tested in a combination (technically no repetition).
other_cols <- c("A", "B", "C", "D", "E", "F")

# Need the same one here for any of the single ones you want to try as random effects.
other_colsRnd <- c("(1|A)", "(1|B)", "(1|C)", "(1|D)", "(1|E)", "(1|F)")

Obviously, there needs to be no repetition of a letter as both a random and a fixed effect.

I have been stuck on this as an algorithm, let alone the coding, for much longer than I had hoped, so any help is much appreciated.

As an example of what I have tried was a variant of this post using dredge:
dredge used on glmer

lmerwrap <- function(formula) { 
    cl <- origCall <-
    cl[[1L]] <-"lmer") # replace 'lmerwrap' with 'lmer'
    # replace "re" with "|" in the formula:
    f <- as.formula("substitute", list(formula, list(re ="|")))))
    environment(f) <- environment(formula)
    cl$formula <- f
    x <- eval.parent(cl) # evaluate modified call
    # store original call and formula in the result:
    [email protected]  <- origCall
    attr([email protected], "formula") <- formula
formals(lmerwrap) <- formals(lme4::lmer)

options(na.action = "") # include at end.  options(na.action = "na.omit")

fm <- lmerwrap(bnParam1 ~ taxon + re(1, D) +  re(1, E) + re(1, F), data = df)
dredge(fm, evaluate = FALSE)

However it does not seem to like to play with lmer in the same fashion, and it would seem impossible (at least for me) to get it to work with the nesting.
I also tried a combo within the combo, but that was such a complete mess, it’s embarrassing.