Edit formula, or convert expanded version into compact notation

Viewed 21

I'm writing a function in R to simplify random effects structures for (generalized) linear mixed effects models that fail to converge (using the lme4 library). I've managed to extract the formula from the model, identify the random effect associated with the least variance, and remove it from the formula. Then, I can join everything back together and use eval(parse(text=formula.string)) to get a formula object to return.

While the formula should work, it's not the most compact notation. For instance, suppose I fit the original model with the following formula:

y ~ x * z + (1 + x * z || group1) + (1 + x * z || group2) 

Then I am able to find that the random effect associated with the least variance is, say, the main effect of z for group2. I successfully construct a new formula that removes only this effect, but then I end up with the following representation (luckily, the dependent variable and fixed effects can be extracted without expansion using lme4's nobars()):

y ~ x * z + 
    (1 | group1) + (0 + x | group1) + (0 + z | group1) + (0 + x:z | group1) +
    (1 | group2) + (0 + x | group2) + (0 + x:z | group2)

This will work, but it's not quite as nice when displaying logging messages regarding what was changed as the ideally more compact notation:

y ~ x * z + (1 + x * z || group1) + (1 + x + x:z || group2)

Is there an easy way to go from the former (either as a formula object or as a string) to the latter without rolling my own function to do it? I could try, but given the high chance I'd make some error somewhere it'd probably not be worth it given that the verbose formula should still work correctly.

Edit: In response to Roland's comment, here's how I'm identifying the random effect associated with the least variance and removing it; I'd be happy if there's a better way to do this, too!

# in context, we only call this if the model has failed to converge
# or if "boundary (singular) fit" occurs,
# based on an any(grepl()) of model@optinfo$conv$lme4$messages
simplify.lmem.formula <- function(model) {
    formula <- formula(model)

    fixed.effects <- deparse(nobars(formula))
    random.effects <- sapply(findbars(formula), deparse)
    
    sds <- as.data.frame(VarCorr(model))
    sds <- sds[order(sds$vcov),]
    least.variance.group <- gsub('\\.[1-9]*$', '', sds[1, 'grp'])
    least.variance.effect <- gsub('\\(Intercept\\)', '1', sds[1,'var1'])
    
    to.remove <- paste0(least.variance.effect, ' \\| ', least.variance.group)
    
    # check to make sure we remove 1 and only 1 effect
    pre.removal.length <- length(random.effects)
    random.effects <- random.effects[!grepl(paste0('(^| )', to.remove, '$'), random.effects)]
    post.removal.length <- length(random.effects)
    if (post.removal.length != pre.removal.length - 1) {
        cat('Error: >1 random effect removed!\n')
        return ()
    }
    
    new.formula.str <- paste0(
                         fixed.effects,
                         ' + (',
                         paste0(random.effects, collapse = ') + ('),
                         ')'
                       )
     # changed from eval(parse(text=new.formula.str))
     # based on Roland's comment
     new.formula <- as.formula(new.formula.str)
}
1 Answers

Try gsub.

fo <- y ~ x * z + (1 + x * z || group1) + (1 + x * z || group2) 

as.formula(paste(gsub('z || group2', 'x:z || group2', fo, fixed=TRUE)[c(2, 1, 3)], collapse=' '))
# y ~ x * z + (1 + x * z || group1) + (1 + x * x:z || group2)
Related