Maybe this is possible already and I'm just missing it, but my impression is that that glm wrappers all basically assume a main terms model and you are using calls glm.fit to fit the models. However, in many applications, we'd like to be able to try out a bunch of different parametric regressions (e.g., with different interactions). As far as I can tell there's not a way to play with the formula interface of the glm function to enable this. Would be appreciated if such a feature could be added (or if I've missed it, please point me in the right direction!).